Weight Function Method With Segment-Wise Polynomial Interpolation to Calculate Stress Intensity Factors for Complicated Stress Distributions

When a crack is detected in components of nuclear power plants during in-service inspections, structural integrity assessment that considers crack growth and failure evaluation has to be conducted in many situations. Integrity assessment procedures have been provided in the ASME Boiler and Pressure Vessel Code Section XI (ASME Code Section XI) [1]. The stress intensity factor (K) is a very important mechanical parameter, and it is necessary to calculate K with high accuracy in order to assess the structural integrity of cracked components with high reliability. Many solutions of K have been proposed and included in codes such as ASME Section XI [1], API 579 [2], RSE-M [3] or JSME (The Japan Society of Mechanical Engineers) Code S NA-1-2008 [4]. Most of these proposed solutions take the stress distribution into account as a polynomial equation that extends to the third or fourth-order. However, many actual stress distributions may be very complicated, such as residual stress distributions or thermal stress distributions in dissimilar materials. In many situations, it is difficult to represent these complicated stress distributions as third or fourth-order polynomial equations over the stress distribution area. Therefore, it is very important to further improve the calculation accuracy of K considering complicated stress distributions.

To calculate a solution of K considering complicated stress distributions, a method based on the influence function method was proposed by Shiratori et al. [5–7]. In this method, the crack surface is divided into several finite elements, and a matrix of influence coefficients is obtained for a specific crack dimension by applying a unit force at every node of finite elements at the crack surface and conducting finite element analyses (FEAs). The results of a large number of FEAs considering different crack dimensions give a series of matrixes of influence coefficients, and these are stored as a database. Although a database with a large number of influence coefficients is required in this method, a value of K considering any arbitrary stress distribution can be calculated by integrating the product of the stress distribution and the influence coefficient.

Recently, a method based on the universal weight function method with piece-wise linear interpolation was proposed by Xu et al. [8] for implementation into ASME Section XI. By using this method, the stress distribution obtained from FEAs or experiments does not need to be curve fitted to a polynomial equation. K is calculated by assuming that the stress varies linearly between two adjacent discrete stress data and utilizing the universal weight function method. Moftakhar and Glinka [9], Anderson and Glinka [10] also investigated discrete stress data and utilizing the universal weight function method. Moftakhar and Glinka represented the discrete stress data as piece-wise linear variation between discrete points. Anderson and Glinka investigated three variations of stress between discrete points: piece-wise constant stress over intervals; piece-wise linear variation over intervals; and piece-wise quadratic variation over intervals. There are many attractive points in these methods. For example, due to the stress distribution is represented as constant, linear or quadratic variation, the exact analytical solution for the weight function integral can be obtained, and then the efficiency of the calculation of K is improved significantly without paying any special attention to the singular point of the weight function and using the complicated numerical integration method for singular points [11]. Also, these methods are easy to apply, and that they do not introduce the uncertainty in the solution of K that results from polynomial approximation of the stress distribution. However, these methods also have some limitations, for example, attention has to be paid to the stress discontinuity [10] and the calculation accuracy for stress distributions with complicated configuration need to be examined.

In this study, a weight function method with segment-wise polynomial interpolation is proposed to calculate K for complicated stress distributions. In this method, the stress distribution area is firstly divided into several segments in accordance with the features of the stress distribution, such as any abrupt change due to the presence of dissimilar materials, and the stress distributions in all segments are curve fitted to segment-wise polynomial equations. Next, K is calculated based on the weight function method and the fitted stress distribution in every segment. Some example solutions for both infinite length cracks and semi-elliptical cracks are calculated and compared with FEA results.

where a is the crack depth, x is the distance from the crack surface, h(x, a) is the weight function, and σ(x) is the stress distribution at the location of the crack.

As shown in Eq. (1), by using this method, a solution of K can be calculated considering independently the stress distribution and the weight function. The weight function depends on the crack shape and geometry of the cracked component. Once the weight function is determined for a certain crack shape in a specified component, K can be calculated for any arbitrary stress distributions.

The development of the universal weight function [14,15] has further increased the usefulness of the weight function method in engineering applications. The universal weight function can represent a wide variety of cracked geometries by using one general form of weight function. For example, for the semi-elliptical crack shown in Fig. 1, the universal weight function for the deepest point (Point A) of the surface crack is expressed as follows:

where M₁, M₂,…, M_n are the coefficients of the universal weight function. The coefficients can be determined for specified crack and component.

Solutions of K based on the universal weight function with four terms have been developed for many crack shapes and structural geometries, such as infinite length or semi-elliptical cracks in cylinders or plates, assuming the stress distribution at the crack location can be represented as a third or fourth-order polynomial equation. Some of these solutions have been included in API 579 [2].

In this study, in order to calculate K considering complicated stress distributions which are difficult to be represented as a third or fourth-order polynomial equation, a weight function method with segment-wise polynomial interpolation is proposed. The concept of the method is illustrated in Fig. 2, where x_j and x_j+1 are the coordinates at the start point and end point of the jth segment.

As shown in Fig. 2, for a complicated stress distribution at the crack location, curve fitting with a polynomial equation over the stress distribution area may not give adequate accuracy. However, if we divide the stress distribution area into several segments considering the feature of the stress distribution, as described later, and curve fit the stress distribution to polynomial equation in every segment, it is easy for the fitted results to give satisfactory accuracy. Therefore, in the weight function method with segment-wise polynomial interpolation, we calculate the solution of K in two steps. The first step is to divide the stress distribution area into several segments and curve fit the stress distribution in every segment to polynomial equation. The second step is to calculate the solution of K based on the universal weight function and fitted results for the stress distributions in each segment.

To date, a universal weight function with four terms has widely been used. In general, the use of more terms is expected to improve the calculation accuracy. In this study, a weight function with six terms was used. The calculation procedure discussed in this paper is still applicable if fewer or more terms of the weight function is chosen.

where a is the crack depth, x is the distance from the crack surface, M₁ to M₅ are the unknown coefficients of the universal weight function.

where N₁ to N₅ are the unknown coefficients.

where $l$ is the length of the crack as shown in Fig. 1.

where G_0B to G_3B are the dimensionless influence coefficients of K at the surface point of the crack.

It can be seen from Eqs. (7a) and (7b) that if the dimensionless influence coefficients G_0A to G_3A and G_0B to G_3B for a specified crack are given, the universal weight function can be determined. To date, many solutions of K have been provided in fitness-for-service codes and fracture mechanics handbooks [1–4,18], for third-order polynomial stress distribution and for important crack shapes and structural geometries. Therefore, using the influence coefficients G_0A to G_3A and G_0B to G_3B as the basis of the universal weight function will not become an obstacle in engineering applications using this method.

where t is the wall thickness of the component, x is the coordinate in the thickness direction, A_ij are the coefficients of ith-order polynomial stress distribution obtained from curve fitting for segment j and n_j is the order of polynomial distribution at jth segment.

If the number of segments is sufficient, the stress distribution with any complicated configuration can be curve fitted appropriately. The purpose of segment division is to ensure that curve fitting gives sufficient accuracy. The candidates for segment division points can be selected so as to correspond to points where there is a trend change (black points) or points where there is an abrupt change (red points), as shown in Fig. 3.

where m is the segment number, and x_j and x_j+1 are the coordinates at the start point and end point of the jth segment, respectively.

The exact analytical solutions of F_ijA and F_ijB can be obtained explicitly in closed-forms from Eqs. (12a) and (12b), by substituting Eq. (3a) into Eq. (12a) and Eq. (3b) into Eq. (12b), respectively. No numerical integration is needed in calculation of Eqs. (12a) and (12b). Therefore, no special treatment such as the complicated numerical singular integration method proposed by Mawarari and Nelson [11] is needed for the singular point for Eq. (3a) at x = a, and for Eq. (3b) at x = 0. Using Eqs. (11a) and (11b), the solution of K can be obtained by representing the stress distributions in different separated segments.

In order to confirm the accuracy of the proposed method, some example calculations for a 360-deg circumferential crack in a cylindrical component were conducted and the solutions of K were compared to the results from FEA.

The analysis model is illustrated in Fig. 4. A cylinder with a circumferential 360 deg. crack at the inner surface was used in the calculation. The wall thickness is 30 mm and the ratio of inner radius to wall thickness R_i/t is 60. abaqus Standard 6.9-1 [19] of Dassault Systems Corp. was employed in FEA. In the analyses, eight-node quadratic axisymmetrical elements were used. In order to represent the singularity of the stress distribution at the crack tip, one edge of the eight-node element connected to the crack tip was collapsed and the midside node was moved to the 1/4 point nearest the crack tip. The example mesh near the crack surface and the mesh near the crack front are shown in Figs. 5(a) and 5(b) for a crack depth with a/t = 0.6. Approximately 2100 elements and 6500 nodes were used in the analyses. The stress distribution considered in the analyses was represented as distribution load acting on the crack surface. Based on the FEA results, the J-integral was obtained at the crack tip. The solution of K was calculated from the results of the J-integral; the third to fifth J-integral paths were used. Plane strain condition was assumed when K was calculated from J-integrals.

For the analysis model shown in Fig. 4, the solutions of K have been provided in C.5.8 of API 579 considering constant, linear, quadratic and cubic stress distributions [2]. In order to confirm the mesh convergence, solutions of K were obtained for a/t = 0.2, 0.4, 0.6, and 0.8, using the FEA model described above. These solutions were compared with the solutions provided in API 579. It was found the solutions obtained from FEA are in good agreement with the provided solutions. The relative difference between them is less than 2%.

Three cases of residual stresses examined in this paper are shown in Fig. 6. All are obtained from FEA [20–22]. The residual stress distribution for case A can be represented as a sixth-order polynomial equation in units of MPa as follows:

On the other hand, the residual stresses for cases B and C cannot be represented with sufficient accuracy over the wall thickness, even if a sixth-order polynomial equation is used, as shown in Fig. 6.

In the example calculations, for the simplicity, only the residual stress was taken into account, although in actual applications there are other stress components such as those due to internal pressure or mechanical bending moment in addition to the residual stress. Therefore, some negative values of K were plotted in the figures as they were obtained.

The weight function method with segment-wise polynomial interpolation was used to obtain the solutions of K for the residual stresses shown in Fig. 6. Based on the principle of segment division shown in Fig. 3, the segment divisions for these residual stresses are shown in Fig. 7. Fourth-order polynomial equations were used for all segments in the three residual stress cases. In case A, the wall thickness was divided into three segments, while the stress distribution area with 0.0 < x/t < 0.8 was divided into three and five segments for cases B and C, respectively. The fitted results for the stress distributions in each segment are shown in Fig. 7 as colored curves. It should be noted that although the stress distribution area for case A can be represented using one segment with a sixth-order polynomial equation, three segments with fourth-order polynomial equations were used. An investigation of different segment division and curve fitting approaches will be provided later, in Fig. 9.

To calculate the solutions of K for complicated stress distributions using the proposed method, it is necessary to utilize the provided influence coefficients G_0A to G_3A, shown in Eqs. (7a) and (7b), as the basis of the universal weight function. The dimensionless influence coefficients used in these calculations are given in Table 1. In the table, the values of the influence coefficients G_0A to G_3A for a/t = 0.2, 0.4, 0.6, and 0.8 were obtained using the FEA model described in Fig. 5. For the dimensionless influence coefficients at a/t = 0.0, the values provided in C.5.8 of API 579 [2] were used.

The results of K using the three residual stress cases shown in Fig. 6 and the influence coefficients given in Table 1 are illustrated in Fig. 8. The solutions using the weight function method with segment-wise polynomial interpolation are shown as solid lines. The FEA results for a/t = 0.2, 0.4, 0.6, and 0.8 are shown as solid circles. For comparison, the solutions of K were also calculated using the conventional method, as expressed by the following equation, representing the residual stress distribution over the wall thickness as a fourth-order polynomial equation:

where σ₀ to σ₄ are the coefficients of the fourth-order polynomial stress distribution, G_0A to G_4A are the influence coefficients. The solutions of K using Eq. (14) are shown as broken lines in Fig. 8.

As previously described, many solutions of K have been provided which consider third or fourth-order polynomial stress distributions. However, many complicated stress distributions, such as those shown in Fig. 6, are difficult to be represented in this manner. If these stress distributions are inappropriately represented as fourth-order polynomial equations when calculating the solutions of K, significant errors will result, as shown by the broken lines in Fig. 8. The maximum absolute differences between the solutions obtained from fourth-order polynomial stress distributions and the solutions obtained from the weight function method with segment-wise polynomial interpolation are 6.0, 5.5, and 11.2 MPa√m at a/t = 0.775, 0.225, and 0.8, respectively, for the three stress cases shown in Fig. 6. The corresponding relative differences are 117%, 20%, and 9%, respectively. Considering the crack growth rate for SCC (Stress Corrosion Cracking) is usually expressed as a power law and the exponent of K, for example, provided in JSME Code S NA-1-2008 [4] is 2.161 for SCC of austenitic stainless steel and 5.186 for SCC of nickel based alloy in a BWR (Boiling Water Reactor) environment, these calculation error in K will significantly influence the predicted crack growth rate and structural integrity assessment.

In contrast, as can be seen from Fig. 8, the solutions using the weight function method with segment-wise polynomial interpolation are in good agreement with the FEA solutions for all three residual stress cases. The relative differences between the solutions obtained from the proposed method and those from FEA for a/t = 0.2, 0.4, 0.6, and 0.8 are less than 0.5% for the three stress cases. These results show that, using only a small number of segments, the proposed method can be applied with high accuracy to the consideration of actual complicated residual stress distributions.

As shown in Eqs. (11a) and (11b), when the universal weight function with segment-wise polynomial interpolation is applied to a complicated stress distribution, there are two approaches that may be taken to segment division and curve fitting in order to improve the accuracy of calculation of K. One of them is to divide the stress distribution area into several segments and use a relatively low order polynomial equation. Another is to increase the order of the polynomial in order to represent the stress distribution with a small number of segments. In order to investigate the difference between these two approaches and provide recommendation to determine how many segments and what is the appropriate order of the polynomial should be used in engineering applications, both approaches were used to calculate solutions of K for the residual stress in case A. Using the first approach, the wall thickness was divided into three segments and curves fitted to the stress distribution in each segment using fourth-order polynomial equations, as shown in Fig. 7(a); in the second approach, a sixth-order polynomial equation curve was fitted to the stress distribution represented as a single segment, as described previously. It can be found from Fig. 7(a) that both approaches gave a sufficiently accurate representation of the stress distribution.

The results for K are shown in Fig. 9. As can be seen from this figure, the solutions of K given by both approaches to segment division and order of the polynomial from curve fitting are in good agreement. This means if the stress distribution can be represented accurately, both approaches to segment division and curve fitting are able to provide sufficiently accurate solutions of K.

Example calculations for a semi-elliptical surface crack in a plate were also performed in order to confirm the accuracy of the proposed method, and the solutions of K were again compared with FEA results.

The analysis model for a semi-elliptical surface crack in a plate is illustrated in Fig. 10. The wall thickness is 30 mm and the width and height of the plate is W = h = 50 t. K was also solved by FEA, 20-node quadratic hexahedron elements being used. Again, abaqus Standard 6.9-1 of Dassault Systems Corp. was employed in FEA. In the analyses, 20-node quadratic hexahedron elements were used. In order to represent the singularity of the stress distribution at the crack tip, singular element shown in Fig. 11 was used at the crack tip. The example mesh near the crack surface is shown in Fig. 12 for a crack depth with a/t = 0.8 and crack aspect ratio with a/$l$ = 0.5. In consideration of the symmetry of the structure, one quarter of the plate was analyzed using around 120,000 elements and 500,000 nodes. Based on the FEA results, the J-integral was obtained for each node along the crack tip. The solution of K was calculated from the results of the J-integral; the third to fifth J-integral paths were used.

For the analysis model shown in Fig. 10, the solutions of K have been provided in RSE-M considering constant, linear, quadratic and cubic stress distributions [3]. In order to confirm the mesh convergence, solutions of K were obtained at crack deepest point for crack aspect ratio with a/$l$ = 0.5, 0.125 and crack depth with a/t = 0.2, 0.4, 0.6, and 0.8, using the FEA model described above. These solutions were compared with the solutions provided in RSE-M. It was found the solutions obtained from FEA are in good agreement with the provided solutions. The relative difference between them is less than 2% [23].

The residual stresses shown in Fig. 6 were again used in these calculations and the segment divisions illustrated in Fig. 7 were also adopted for the semi-elliptical surface crack, with fourth-order polynomial equations being employed for all segments in the three residual stress cases. For the analysis model shown in Fig. 10, the residual stress was considered to vary in the wall thickness direction only. From the viewpoint of simplicity, the residual stress was assumed to be constant in the width direction of the plate. The stress distribution was represented as distribution load acting on the crack surface in the analyses.

As the basis of the universal weight function used in the calculations, the dimensionless influence coefficients at the deepest point of the crack are given in Table 2 for a crack with aspect ratio a/$l$ = 0.5 and in Table 3 for a crack with aspect ratio a/$l$ = 0.125.

In Tables 2 and 3, the values of the influence coefficients G_0A to G_3A at a/t = 0.2, 0.4, 0.6, and 0.8 were obtained using the FEA model illustrated in Fig. 12. For the dimensionless influence coefficients at a/t = 0.0, the values provided in RSE-M [3] were used.

The results of K at the deepest point of the semi-elliptical surface crack for the three residual stress cases illustrated in Fig. 6 and the influence coefficients given in Tables 2 and 3 are shown in Fig. 13. The solutions using the weight function method with segment-wise polynomial interpolation are represented as solid lines. The FEA results for a/t = 0.2, 0.4, 0.6, and 0.8 are shown as solid circles. For comparison, the solutions using the following equation by representing the residual stress distribution as fourth-order polynomial equation over the wall thickness are shown as broken lines:

As with the solutions for a 360 deg. circumferential crack at the inner surface of a cylinder, the solutions of K obtained by expressing the residual stress distribution as a fourth-order polynomial equation over the wall thickness are of insufficient accuracy when compared with the FEA solutions. The maximum absolute differences and the corresponding relative differences between the solutions obtained from fourth-order polynomial stress distributions and the solutions obtained from the weight function method with segment-wise polynomial interpolation are summarized in Table 4 for the three stress cases shown in Fig. 6. The reason for this significant difference is that the complicated stress distributions shown in Fig. 6 cannot be expressed with adequate accuracy as fourth-order polynomial equations.

In contrast, as can be seen from Fig. 13, the solutions using the weight function method with segment-wise polynomial interpolation are in good agreement with the FEA solutions for all crack dimensions in all cases of residual stress. The relative differences between the solutions obtained from the proposed method and the FEA solutions are less than 0.5% for all crack dimensions in all residual stress cases. These results show that the proposed method can also be applied with high accuracy to semi-elliptical surface cracks giving consideration to actual complicated residual stress distributions.

In this study, a weight function method with segment-wise polynomial interpolation was proposed for the calculation of the stress intensity factors for complicated stress distributions. Some example solutions for both infinite length cracks and semi-elliptical surface cracks were compared with the corresponding FEA results. Based on the numerical investigation, it can be concluded that the stress intensity factor can be obtained with high accuracy using the proposed method, the actual complicated stress distributions being successfully taken into consideration for both infinite length cracks and semi-elliptical surface cracks. If the stress distribution can be expressed accurately, different approaches to segment division and curve fitting can be adopted without influencing the accuracy of the solutions. For general complex stress distributions, it is not necessary to employ a large number of segments.