Interaction of surface cracks subjected to non-uniform distributions of stress

Closely-spaced surface cracks in structures interact with each other when subjected to load. The degree of interaction depends strongly on the distribution of stress that is applied. In pressure boundary components, thermal shock, residual stress and global bending can all cause load distributions that are non-uniform through the wall thickness. A wide range of crack pairs subject to various non-uniform stress distributions have been modelled using finite element analysis. Cracks sometimes interact more strongly under non-uniform loading than when loaded in uniform tension. Consequently, interaction criteria developed by considering uniform tension may not be inherently conservative for all loading conditions. A simple weight function method for determining the interaction of twin cracks under an arbitrary through-wall stress is presented, and weight function coefficients for a wide range of crack sizes and aspect ratios are given. © 2017 The Author. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Most procedures for structural integrity assessment include rules for establishing whether two or more closely-spaced defects in a structure interact. When performing an integrity assessment, it is often useful to know whether failure mechanisms that are initiated by the presence of a defect are influenced by the presence of other defects nearby. This tells the assessor whether the defects can be considered individually or whether they will have to be considered jointly. If the interaction is significant and the defects need to be jointly considered, this can increase the complexity of the analysis. It can also introduce additional conservatism if the interaction can only be accounted for in an approximate manner.
Many structural integrity assessment procedures, including R6 [1] and BS 7910 [2], include simplified interaction criteria that are largely based on Linear Elastic Fracture Mechanics (LEFM). It is assumed that if two crack-like defects are close enough for the stress intensity factor at each to be affected significantly by the presence of the other, then other failure mechanisms may be affected likewise. For example, fatigue crack propagation under Mode I loading is normally assumed to be controlled by the variation in K I . Therefore, the fatigue crack growth rate at one crack can be affected by the proximity of another and this can affect how the cracks coalesce [3e5].
Two-dimensional cracks in the through-thickness plane (as shown in Fig. 1) are a particularly important concern for the integrity of pressure vessels. The orientation of such defects means that they can be subjected to opening-mode stresses due to pressurisation, thermal shock and welding/cladding residual stresses. Furthermore, defects of this type can be created by solidification cracking in welds during manufacture or by environmentallyassisted cracking during service. These cracks can be subjected to distributions of stress which vary significantly through the thickness of the pressure vessel wall (Fig. 1). Stress distributions that vary through the wall thickness can be caused by thermal shock, residual stress or local bulging of the pressure vessel. In thickwalled pressure vessels the circumferential stress due to pressurisation is inherently non-uniform, and is greatest at the internal surface.
Although it is known that adjacent cracks in a linear-elastic material interact to a greater or lesser degree depending on the through-thickness stress distribution, all widely-used code criteria for judging whether two co-planar cracks will interact are based purely on the geometry of the defects and are independent of their loading mode. This includes the criteria used by the ASME B&PVC Section XI [6], JSME S NA1 [7], API 579e1 [8], R6 [1], BS 7910 [2], GB/ T 19624e2004 [9] and SSM 2008:01 [10] assessment procedures. The effect of through-wall variations in stress on crack interaction has been investigated by several groups of researchers. Using the body force method, Murakami and Nemat-Nasser [11] first showed that the level of interaction between two coplanar semi-elliptical surface cracks was significantly different under pure bendingtype loading compared with uniform tension. More recently, a pair of studies by Sethuraman et al. [12] [13], calculated interaction factors for pairs of twin semi-elliptical surface cracks in finitethickness plates, again showing that the level of interaction could be different. Finally, Coules [14] performed calculations for a large number of pairs of dissimilar surface cracks, concluding that the effect of through-thickness stress variation in increasing interaction could cause defects to be wrongly classified as non-interacting by commonly-used interaction criteria.
This article examines the effect of through-wall variation in stress on the elastic interaction of crack-like defects. A method for estimating the interaction factor based on the through-wall stress for pairs of semi-elliptical surface cracks is presented.

Method of weight functions
For a single semi-elliptical surface crack subjected to a onedimensional distribution of stress, the stress intensity factor as a function of position on the crack tip line can be expressed as: where K I is the Mode I stress intensity factor, f is a parametric angle defining the position on the crack tip line (see Fig. 2), a is the crack depth and the superscript d ¼ ∞ denotes a crack located remote from any other. The weight function mðy; fÞ relates the position y of a pair of unit normal line loads to the resulting stress intensity factor K I ðfÞ, and it depends only on the geometry of the crack and the solid body containing it. The crack itself can be defined in terms of the crack depth a, crack width 2c and wall thickness t as shown in Fig. 2. It is assumed that the applied normal stress s zz varies in the through-wall dimension (y) only. The ratio of the stress intensity factors for an interacting crack and a remotely-located (but otherwise identical) crack defines the interaction factor g: where the superscript int denotes crack that is in proximity to another. Hence: Since the weight functions m int and m d¼∞ are always nonconstant functions of y [15], the interaction factor always depends on the through-wall distribution of stress s zz ðyÞ.

Deepest point
To remove the need for explicit integration over the crack length in Equation (1), some schemes for estimating K I approximate the through-wall distribution of stress s zz ðyÞ using a polynomial [1], [16]. Here, the distribution of stress over the wall thickness t is approximated using: where a is the crack depth, M 1…3 are geometric constants, and the subscript A indicates coefficients and functions associated with the  deepest point on the crack tip line (see Fig. 2b). This functional form can also be used to approximate the weight function for the deepest point in an interacting pair of twin semi-elliptical surface cracks ( Fig. 2): The constants M 1…3 for a given crack geometry can be determined using two reference values of K I for different through-wall loading conditions [18], [19]. Equation (2) can now be solved to give g in terms of s i : where the coefficients A ij are given in Table 1 (7), the interaction factor can be found straightforwardly from the non-uniform stress distribution s zz ðyÞ without explicit integration over a. Determining the coefficients M d¼∞ 1A , M d¼∞ 2A and M d¼∞ 3A requires two numerical models of a single surface defect: one in which the crack is subjected to uniform tension and one where it is subjected to a linearly-decreasing stress distribution with a finite stress at the surface and zero stress at the deepest point [18]. Similarly, determining the coefficients M int 1A , M int 2A and M int 3A requires numerical models of the twin surface cracks under the same two loading conditions.

Surface point
For points close to the intersection between the tip of a single semi-elliptical defect and the surface, Shen & Glinka [18] proposed a weight function in the following form: where the subscript B indicates values/functions associated with a surface point. As above, the same functional form can be used for a surface point in a pair of twin semi-elliptical surface defects. The weight function for a surface intersection point in-between the two cracks is: Equation (3) can now be solved in the same manner as Equation (6) to give: for the surface point. The coefficients B ij are given in Table 2 and For materials with ns0, at the intersection between the crack tip line and the free surface the strength of the stress singularity deviates from s ij ðqÞfr À 1 2 [18], [20]. This creates a boundary layer close to the free surface where the crack tip stress field is not characterised by K I . However, while the use of a weight function for K I to characterise the crack tip stress field at this location is theoretically unsound, the ratio g can still be used to indicate the general effect of crack proximity on the stress field which does occur.

Model description
The finite element method was used to calculate the Mode I stress intensity factor as a function of position on the crack tip for twin pairs of proximate defects and for single defects with the same dimensions. Models were created using the Abaqus/CAE v6.12 finite element pre-processor [21] and solved using Abaqus/Standard v6.12 [22]. To generate the large number of parametric models required for this study, the Abaqus/CAE Python scripting interface was used in conjunction with custom code written in MATLAB [23]. This process is described in more detail in a previous paper [14].
For pressure vessels where the wall thickness is much smaller than the vessel diameter, the wall can be geometrically approximated using a flat plate. A plate containing the cracks was modelled with unit thickness, a width of 100 units, and a total length of 200 units in the loading direction (see Fig. 3). Due to the reflective symmetry of the plate about the plane containing the cracks, only one half of it was modelled. The regions surrounding the crack tips were meshed using 8-noded linear brick elements, with 6-noded linear wedge elements adjacent to the crack tip line itself. Surrounding this was a region containing 10-noded quadratic tetrahedral elements, while the remainder of the plate was meshed using 8-noded linear brick elements. This meshing arrangement is shown in Fig. 3. All cracks had 81 equally-spaced nodes on the crack tip line. The models were loaded using an equivalent non-uniform pressure applied to the faces on the cracks [14]. The plate material The strain energy release rate was extracted from the finite element results using the equivalent domain integral method [24]. Since the cracks are under Mode I loading only, the strain energy release rate can be related directly to the Mode I stress intensity factor. For the intersection between the crack tip line and the plate surface, the crack tip node one element away the free surface was used for stress intensity factor extraction. Using results from both the single crack and the twin cracks, the interaction factors for the deepest and surface-intersection points were calculated according to Equation (2).
The coefficients needed for the weight function method described in Section 2.2 were determined from stress intensity factor results for cracks subjected to uniform tension and pure bending. First, superposition was used to determine the stress intensity factor for a linearly decreasing stress distribution over the crack length from the results for tension and pure bending: K L:D:

Stress distribution test cases
To test the accuracy of the weight function method, additional finite element models were run using three 'test-case' stress distributions listed in Table 3. For each of these loading cases, all combinations of crack depth, aspect ratio and separation distance listed in Table 3 were modelled. The interaction factors were calculated directly from the finite element results and then compared with interaction factors calculated via the weight function method described in Section 2.2.
The first stress distribution listed in Table 4 is taken from the structural integrity assessment procedure R6 [1]. In that procedure, it is used as a conservative estimate of the distribution of transverse residual stress in austenitic and ferritic steel plate butt welds. It was constructed as an approximate upper bound to a range of measured and modelled residual stress distributions for such welds. The second test case is the predicted distribution of hoop stress in the wall of a WWER 440 Model 213 reactor pressure vessel, 1200 s from the start of a thermal shock event initiated by inadvertent opening of the pressuriser safety valve. This distribution was calculated as part of an IAEA benchmarking study in which six participating organisations independently predicted the through-wall stress distribution using FEA, with good agreement [25]. The 5 th -degree polynomial used to represent it fits closely to the modelled data. The final stress distribution is a simple combination of tension and reverse bending, such that the stress always increases between the surface y ¼ 0 and the crack tip.

Interaction factors
An example result from finite element analysis of the semielliptical cracks is shown in Fig. 5. Fig. 5a shows the Mode I stress intensity factor as a function of position on the crack tip line for a single crack with dimensions a t ¼ 0:5 and a c ¼ 1. Fig. 5b shows the result for a pair of such cracks with a separation distance of For twin cracks of a given crack size and aspect ratio in uniform tension, the interaction factors at both the deepest point and surface point always increase with decreasing crack separation. Fig. 6 shows stress fields around pairs of cracks with decreasing separation distance; the stress field in-between the twin cracks becomes more intense as the cracks approach each other. The interaction factor at the surface point (which is the point of closest approach between crack tip lines) is always greater than the interaction factor at the deepest point. This is illustrated in Fig. 7 and is in agreement with previous studies by Sethuraman et al. [13] and Coules [14].

Weight function method accuracy
The weight function coefficients for single semi-elliptical surface cracks (M d¼∞ iA , M d¼∞ iB ) and twin pairs of such cracks (M int iA , M int iB ) were calculated from the finite element results. Coefficients for the full range of crack pair geometries investigated here are provided in the supplementary data. When used in Equation (7) and Equation (10), these coefficients enable calculation of the interaction factor for twin pairs of semi-elliptical surface cracks subjected to arbitrary through-wall distributions of stress. Fig. 8a shows a comparison of interaction factors calculated via FEA with those calculated using the weight function method described in Section 2.2, for the deepest point on pairs of interacting cracks. 140 crack pairs are shown, each loaded with the three stress distributions shown in Table 4 and Fig. 4. There is excellent agreement between the interaction factors calculated using the two different methods. As can be seen from Fig. 8b, the error is generally less than 1%, however up to 5% error was observed for a few geometries where the stress intensity factor at the deepest point was close to zero. Likewise, Fig. 8c and d shows the accuracy of the weight function result for the surface point: the error with respect to the FEA result is less than 1% in all cases.

The effect of stress distribution on stress intensity interaction
One simplifying method that can be used to determine the effect of a non-uniform distribution of stress on the stress intensity factor at a crack is to decompose the stress distribution into tension, bending and self-equilibrating components [26]. This technique is used in the assessment procedure BS 7910, for example [2]. The contributions of these three components to the stress intensity factor are additive: they can be calculated separately and then summed to produce the overall K I result. However, the interaction factor g is a ratio of stress intensity factors that are different functions of s zz ðyÞ (see Equation (3)) so the contribution of any polynomial component of the stress distribution to the interaction factor depends on the other components. So the tension or bending parts of a stress distribution do not necessarily contribute more to the interaction factor than any self-equilibrating component. However, from the form of Equation (7) and Equation (10) it is evident that higher-order polynomial components of the stress distribution will generally contribute less to each part of the interaction factor ratio. This is due to the presence of the factor a t i and the decreasing magnitude of the coefficients A ij and B ij with increasing i (see Tables 1 and 2).
For the interacting semi-elliptical surface defects investigated here, the point of closest approach between the cracks occurs at the intersection between the crack tip line and the surface of the plate (see Fig. 2a). Generally, the interaction factor is greatest at this point of closest approach. For single cracks of low aspect ratio ( a c < 1) subject to a uniform stress distribution, the highest stress intensity factor always occurs at the deepest point on the crack tip line. For pairs of such cracks, a high interaction factor at the surface point has no effect on fracture initiation, which will occur elsewhere. However, for high aspect ratio cracks ( a c ! 1) under uniform tension and for cracks of any aspect ratio subject to bending, residual stress or thermal shock, the stress intensity factor at the surface point is often greater than at the deepest point (see Fig. 5). Therefore, under these conditions the region of highest K I for a single crack coincides Table 4 Polynomial coefficients defining three through-thickness stress distributions (see Equation (4)) that were used for testing the method of interaction factor determination. with the region where there is greatest interaction risk when two cracks occur. Putting aside considerations of constraint loss at the free surface, it can be concluded that the risk of interaction is generally greater when high aspect ratio cracks and/or non-    uniform loading are involved, than it is for low aspect ratio cracks under uniform tension.
An example of this effect is shown in Fig. 9. The results are for a pair of twin cracks with a t ¼ 0:5, a c ¼ 0: _ 3 and d t ¼ 0:5. From its geometry, this crack pair would be classed as having no significant interaction by the current R6/BS 7910 interaction criterion [1], [2]. We can define a 'global' stress intensity factor as [14]: which is the ratio by which the maximum stress intensity factor found anywhere on either crack is increased by the proximity of the cracks to one another. This is illustrated in Fig. 9a. The factor g G can be used as a measure of how significant interaction is in controlling brittle fracture initiation. Although the crack pair in Fig. 9 shows a high level of interaction close to the surface point when under uniform tension, the stress intensity factor is always lower here than at the crack's deepest point. However, with a non-uniform stress distribution the point of greatest interaction coincides with the point of greatest stress intensity factor for a single crack. In this case we have g G ¼ 1:159, i.e. an increase in the maximum stress intensity factor of 15.9%.

Implications for assessment methods
In a recent study [14], it was shown that non-uniformity of the stress distribution can strongly affect the interaction factor, and that this effect can cause unintended non-conservatism in the criteria used for determining the significance of defect interaction in integrity assessment procedures. Interaction criteria based on stress intensity factor results for uniform tension may not be conservative when applied to cracks that are acted on by a nonuniform stress distribution. However, the weight function method described in Section 2.2 provides a straightforward means of determining the interaction of surface defects that are loaded nonuniformly. This could be used as the basis for guidance in dealing with strongly non-uniform through-wall distributions of stress.
For example, if the through-wall distribution of stress applied to the flaws is uniform then it is normally sufficient to use the interaction criteria for coplanar surface flaws detailed in R6 Figure II.3.6 [1] or BS 7910 Figure 12 [2]. Otherwise, if the through-wall distribution of stress is non-uniform then a sensible approach would be to: 1. Re-characterise the defects as a pair of twin cracks, with each crack large enough to fully enclose each defect. The spacing of the cracks (d) should be the same as the observed spacing of the defects. 2. Using Equation (7) and Equation (10), determine interaction factors for the surface and deepest points when the twin cracks subjected to the anticipated (non-uniform) distribution of stress. The through-wall stress distribution should be approximated as a 5 th -degree polynomial.
Then, interaction is judged to be insignificant only when both of the following inequalities are satisfied: 1 g A < g crit 1 g B < g crit (13) where g crit is a critical interaction factor representing the allowable increase in stress intensity factor. A typical value of g crit that has been used for formulating interaction criteria is 1.1, but this has no rigorous basis and lower values have also been used. Although Equation (13) only considers the crack's deepest point and surface intersection point, under most realistic loading cases the greatest stress intensity factor will occur close to one of these locations.

Limitations
All results presented here deal with crack interaction under linear elastic conditions. The interaction of defects subject to other failure processes such as plastic collapse and creep has not been quantified in this study. However, since plasticity and creep deformation tend to reduce the effect of one stress concentration on another, considering elastic interaction only in these cases is generally believed to be conservative [1], [27]. Likewise, additional loss of crack-tip constraint due to the proximity of adjacent defects, and the proximity of any free surfaces, will allow structures to bear higher loads to failure than is implied by the single-parameter approach to characterisation of the crack tip field used here. Therefore, this approach gives an inherently conservative estimate Fig. 9. Interaction of a pair of cracks ( a t ¼ 0:5, a c ¼ 0: _ 3, d t ¼ 0:5) subjected to uniform tension and the R6-presecribed plate butt weld transverse residual stress distribution (see Fig. 4). a.) Normalised stress intensity factor, b.) interaction factor. of the level of interaction that will occur for a real elastic-plastic material.
The weight function method presented in Section 2.2 does not consider the possibility of contact between the opposing faces of a crack. Crack face contact may modify the stress intensity factors and hence the interaction factor for a given crack pair. Where the loading mode is strongly non-uniform, with compressive stress acting over part of the prospective crack plane, the possibility of crack face contact is greater than for cracks under uniform tension. When crack face contact is not considered, negative stress intensity factors can be calculated over part of the crack tip line leading to unrealistic interaction factor results. An example of a pair of deep ( a t ¼ 0:75) cracks is shown in Fig. 10. When a single deep crack or a pair of such cracks is subjected to through-wall bending, a negative stress intensity factor is calculated for the deepest section of the crack tip line (see Fig. 10a). This is not physically meaningful and results in an interaction factor distribution containing asymptotic discontinuities as shown in Fig. 10b. To account for this effect a lower bound on g has been included in Equation (13), ensuring that any combination of crack pair and loading that shows this effect cannot be assumed to be non-interacting.

Conclusions
1. Elastic interaction between coplanar surface cracks is strongly affected by the through-thickness distribution of stress applied to them. 2. It is not inherently conservative to assume that closely-spaced cracks subjected to bending, residual stress or thermal shock will interact to the same degree as they would under uniform tension. 3. A weight function method can be used to calculate interaction factors for non-uniform loading accurately. Weight function coefficients have been determined for pairs of identical semielliptical surface cracks subjected to any arbitrary throughwall distribution of stress that can be closely approximated with a 5 th -degree polynomial. 4. Several structural integrity assessment procedures contain criteria for judging the significance of defect interaction. These criteria are normally independent of the loading mode. Care should be taken when applying such criteria to cracks that are subjected to strongly non-uniform loading.