A geometry-dependent generalized shape function for calculation of stress intensity factor for axially cracked thin-walled tubes

Geometric shape function, f(a/W), is a correction factor that accounts for finite specimen size in calculation of stress intensity factor (SIF). Pin-loading tension method is suitable for fracture toughness testing of axially cracked thin-walled tubular specimens. But unique expression for f(a/W) of such specimens does not exist as it is sensitive to tube geometry. In this work, a generalized form of geometric shape function has been derived through a detailed finite element analysis and the outcome is validated with actual experimental results. Using the general shape function, SIF for any axially cracked thin-walled tube having dimensions within the range can be predicted.


Introduction
Thin-walled tubes such as fuel cladding tubes made of zirconium alloys, steam generator (SG) tubes made of incoloy-800 etc. find wide applications as structural members in nuclear industry. Cladding tubes encapsulate radioactive fuel pellets and isolate them and other harmful fission products from the surrounding environment by providing physical barrier. During service and storage, they are subjected to neutron irradiation, internal gas pressure due to fission gas release, a highly corrosive environment leading to oxidation and hydrogen degradation, pellet-clad mechanical interaction and high temperature (Meyer et al. 1996;Nikulina 2003;Leclercq et al. 2008). Again, steam generator tubes face flow-assisted vibration, corrosion, and differential thermal stress due to wide difference in temperature between the inner and outer walls (Saint Paul and Slama 1992). Synergistic effect of all these factors turns the thin-walled tubular components more prone to axial cracking due to dominance of circumferential hoop stress. For longer uninterrupted plant operation, residence time of these components is to be maximized within the plant and that necessitates a thorough structural integrity analysis of such tubes. The geometry of the tubes does not allow one to resort to any direct and straightforward standard procedure by which a valid fracture property can easily be evaluated. Out of the various non-standard test procedures, Grigoriev et al. (1995) reported a novel test method called Pin-loading tension test (PLT) in which a section of the thin-walled tube itself can be loaded in combination with suitable jig and fixture assembly ( Fig. 1) to generate meaningful fracture data. Samal et al. (2010a, b) in recent past customized the design of the jigs for cladding tubes used in Indian reactors and evaluated the various geometric factors such as geometric shape function f(a/W) required for calculation of stress intensity factor (SIF), g and c functions required for calculation of J integral from experimental data with a detailed agenda involving experiment and finite element method (FEM). Using those functions, they also calculated J-R curves for two different cladding tubes used in Indian nuclear reactors with experimental data both by load normalization (Sanyal et al. 2011a, b;Samal and Sanyal 2012;Sanyal and Samal 2012a, b) and load separation procedure and compared the results (Samal et al. 2011). In recent past, Sanyal and Samal (2012a) also assessed axial cracking of steam generator tubes used in Indian pressurized heavy-water reactor (PHWR) for generation of electricity with the same method.
Till date, the authors evaluated the geometric shape function, f(a/W), both experimentally and using FEM for seven different thin-walled tubes used in different nuclear reactors. The function f(a/W) is important as it takes care of the effect of finite specimen size in calculation of SIF according to the equation (Samal et al. 2010a, b), Where E is young's modulus, C is elastic compliance for a given crack length, a, W is specimen width (= 19 mm for all current cases) (more specifically, W is width of the specimen-fixture assembly), t is tube wall thickness, and K I is SIF in mode I loading. Here, it is to be noted that for PLT specimens, a/W is essentially the ratio of sum of specimen crack length (a) and virtual crack length of the fixture (b = 8 mm for all current cases as W = 19 mm) to W (Grigoriev et al. 1995;Samal et al. 2010aSamal et al. , b, 2011Sanyal et al. 2011a, b;Samal and Sanyal 2012;Sanyal and Samal 2012a, b). In other words, for calculation of a/W, the crack length due to specimen-fixture assembly has to be considered rather than considering crack length of the specimen alone. Using K I , elastic contribution in calculation of J integral (J el ) is calculated using the following equation for the purpose of construction of J-R curve (Sanyal et al. 2011a, b;Samal and Sanyal 2012;Sanyal and Samal 2012a, b;Samal et al. 2011): It is noticed, since the past experience that the function, f(a/W) is sensitive to tube geometry, i.e., as the combination of tube diameter (D) and wall thickness (t) changes, the function f(a/W) also changes. From the analogy of the loading mechanism between a PLT specimen and a compact tension (CT) specimen, it may apparently be aspired that the functional form of f(a/W) for them should be identical; in other words, it can be thought of that f(a/W) should be independent of geometry for PLT specimens as it is so for CT specimens. But as it requires two geometric parameters, i.e., D and t to define the geometry of a thin-walled tube, the function f(a/W) is found to depend on both D and t, and there is no report available in open literature so far to the best of authors' knowledge to correlate f(a/W) with D and t for PLT specimens. If f(a/W) can mathematically be associated with D and t then one can derive the functional form of f(a/W) for thin-walled tubes with any combination of D and t and also can evaluate the value for f(a/W) for any known a for any thin-walled tube. Thus, establishment of such a function can save much effort needed in experimental and analytical derivation of f(a/W) for axial splitting of any tube having a specific geometry and thus, calculation of SIF and J el would be easier.
In the present work, f(a/W) function is found for tubes with different D and t over a wide range through FE analysis. After generation of data, the trend of variation of f(a/W) with change in D and t is observed and accordingly the f(a/W) data are correlated with D and t using suitable fitting functions with the method of least square adapting a multiple-step regression analysis and a generalized geometric shape function [f(a/W, D/W, t/W)] is suggested. The resulting function is also validated with available experimental results. The function can predict the value of f(a/W) for a given set of geometric parameters within the range, eliminating the necessity of resource-intensive experimental procedure and FE analysis practiced earlier.
This paper is organized in five sections. In the section titled 'Analytical framework', the details of FE modeling and simulation are described followed by the discussion about the method adopted for correlation of f(a/W) with D and t. The results are produced followed by critical discussions in the next two following sections with emphasis in validation of the model with experimental results. Finally, the concluding remarks with future scope of similar works are suggested in a separate section and the necessary input for derivation of the experimental data is presented at the ''Appendix''.

Analytical framework
The FE modeling and simulation The actual picture of the test setup used for carrying out PLT test is shown in Fig. 1. The details of the specimen  (Grigoriev et al. 1995;Samal et al. 2010a, b;Sanyal et al. 2011a, b;Samal and Sanyal 2012;Sanyal and Samal 2012a, b). As shown in Fig. 1, the specimen is a 13-mm long section of a thin-walled tube that essentially contains four coplanar notches (diametrically opposite two notches at both end; the loading end contains sharper notches and the other end contains thicker rectangular slots for avoiding compressive load at that end) and it is coupled with the fixture having two semi-cylindrical mandrels capable of loading the specimen through inner surface in a mode I loading fashion analogous to CT specimens; the axis of mutual rotation of the fixture halves is determined by a small pin inserted within the groves of the halves. As per the specimen and fixture design, the entire PLT test setup (i.e., the tubular specimen coupled with the mandrel) has been discretized in 3D domain for FE analysis. This is required as the load is applied at the free end of the mandrel to load the cracked tubular specimen in mode I fashion. Keeping in mind about the earlier experiments and analyses done, to cover a wide range of tube geometry, thin-walled tubes having diameter (D) varying from 4 mm to 15 mm with an interval of 1 mm and wall thickness (t) varying from 0.3 mm to 1.2 mm with an interval of 0.1 mm are considered in the present study. Hence, altogether 12 9 10 different combinations of D and t are presented in the current case. In each case, seven different notch lengths (a), 1 mm to 7 mm with an interval of 1 mm are considered as has been done earlier (Samal et al. 2010a(Samal et al. , b, 2011Sanyal et al. 2011a, b;Samal and Sanyal 2012;Sanyal and Samal 2012a, b) in individual cases. Hence, altogether 7 9 120 different FE models are generated prior to simulation. A representative model is shown in Fig. 2. The mechanical properties of the material used in the analysis are for a zirconium alloy and are given in Table 1. The mechanical property of the mandrel material (i.e., steel with Young's modulus of 210 GPa and Poisson's ratio of 0.3) is different from that of the tubular specimen.
One-quarter of the test setup has been modeled for the crack and loading configuration and symmetry in the specimen geometry. The symmetric plane of the tubular specimen and the same of the mandrel is fixed in z-direction only for the symmetric boundary condition. A preset displacement is applied through displacement control mode at mandrel end in y-direction and the motion of the loading point is constrained in the x-direction (i.e., axial direction of the tube) as this point is attached to the loading axis of the test setup. The un-cracked ligament of the tubular specimen is fixed in y-direction (but it can move in x-and z-directions) considering the symmetry condition. The line of the loading pin (i.e., the fulcrum) is also fixed in ydirection only. The FE discretization consists of 20-noded iso-parametric brick elements with 3 9 3 9 3 Gauss point integration (full integration scheme) which is suitable in avoiding spurious locking of the FE mesh as deformation proceeds. Bathe (1995) discussed in detail about these elements. It was observed by the authors (Samal et al. 2010a, b) that use of 20-noded elements gave better result than using 27-noded brick elements. Each FE node without any boundary condition is free to move in the x-, y-, and zdirections. For using optimum number of 3D elements in FE mesh, fine discretization has only been used near the crack-tip for all the specimens to reproduce the stress gradient near the crack-tip. The sizes of the 3D elements are 0.2 9 0.2 9 0.2 mm 3 near the crack-front (Fig. 2). For each of the different 120 cases with different combinations of D and t, seven different meshes are generated corresponding to seven different crack lengths, i.e., a varying from 1 to 7 mm by altering the position of the refined mesh along x-direction as shown in Fig. 2 and changing the mirror boundary conditions at the crack-plane.
An in-house FE code is used for making the simulation. For finding the linear elastic SIF (i.e., K I ) for this experimental setup, the data from elastic analysis are essential. For evaluating the SIF, the analysis is carried out for the applied displacement of 0.01 mm in each specimen and the corresponding compliance values are calculated for each value of a/W. Many approaches exist in the literature (Paris and Sih 1965;Ingraffea and Manu 1980;Raju and Newman 1977;De Lorenzi 1982;Parks 1974) to determine the SIF values. If all approaches are broadly categorized, there  are two different methods with which SIF can be derived from the compliance data. The direct method that utilizes the result of the general analytic solutions of crack problems follows either from the stress-field or from the displacement-field extrapolation. The extrapolation techniques are based on asymptotic displacement and stress fields near the crack-front (Paris and Sih 1965): by the displacement-field approach (Ingraffea and Manu 1980) and by nodal-force approach (Raju and Newman 1977). The indirect method on the other hand is one where K is determined via its relation with other quantities such as the compliance, the elastic energy (G), or the J integral (De Lorenzi 1982;Parks 1974). For indirect methods, it is possible to calculate the compliance for a range of different crack sizes; from the data using numerical differentiation with respect to crack size SIF can be evaluated. Alternatively, the total elastic energy (U) contained in the specimen can be found for a number of crack sizes. Numerical differentiation of U with respect to crack size (a) gives the elastic energy release rate (G), from which K can be evaluated as K 2 = EG. Finally, for the elastic case, since G = J, SIF can also be calculated from J integral, derived using domain integral method. The displacement extrapolation method is adopted in this work, where the SIF can be written as (neglecting higher order terms) where K I is SIF in the opening mode of the crack, v is the opening displacement ahead of crack-tip, r is the distance from the crack-tip, E 0 = E/(1 -m 2 ) for plane strain loading condition whereas, E 0 = E for plane stress loading condition, E and m are the elastic modulus and Poisson's ratio for the material, respectively. v is the nodal displacement and it is measured for nodes residing both at inner and outer surface of the specimen (Fig. 2) at the crack-front. Hence, there would be one K I corresponding to inside node and another corresponding to outside node for each case. In subsequent analysis, both cases, i.e., inside and outside are independently dealt followed by another case where the average result of the two is considered. From SIF, the geometric shape function, f(a/W) value for all case is derived using Eq. 1.
Correlation of the geometric parameters with the derived shape function The f(a/W) values for all 120 cases are plotted against a/W individually for each t with all Ds to observe the variation of the function with D and t. A representative plot is shown in Fig. 3 for inside nodes in the case of t = 0.7 mm for all D. In a relative sense, it can be said that f(a/W) is more sensitive to change in diameter rather than that in wall thickness. That can be confirmed if one plots f(a/W) vs. a/ W for a certain D for all t (Fig. 4). Same has also been proved earlier by Chakravartty et al. (2009) in a smaller scale with only three different real thin-walled cladding tubes. From the given data, it can be found that the sensitivity of f(a/W) on D for a given t is more if D is less (i.e., for D = 4 mm, 5 mm etc.). But on increasing D, i.e., for D [9 mm, f(a/W) is relatively less sensitive to change in D. For finding any correlation between f(a/W) and D, all curves plotted for f(a/W) vs. a/W are fitted and a third-order polynomial in a/W in the following form is found to impart best fit in all cases: f ða=WÞ ¼ a 0 þ a 1 ða=WÞ þ a 2 ða=WÞ 2 þ a 3 ða=WÞ 3 for a fixed D and t   In each case, a regression coefficient of C0.999 is achieved. For all the 120 cases pertaining to different combinations of D and t, different values of the coefficients a 0 , a 1 , a 2 and a 3 are found. It can be said that to find a function f(a/W, D) from the data for a given t, the coefficients a 0 , a 1 , a 2 and a 3 must be expressed in terms of D/W (here D is divided by W to impart a non-dimensional form like a/W). To accomplish that, the coefficients a 0 , a 1 , a 2 and a 3 are individually plotted against D/W for each t and the variation of a 0 and a 1 with D/W is shown in Figs. 5 and 6 for all t. Variation of a 2 and a 3 with D/W is similar to that for a 0 and a 1 , respectively. Each of the curves is fitted and a fifth-order polynomial in D/W in the following forms is found to give the best fit for all cases: Equation 5 is valid for a given t. In each case, a regression coefficient of C0.99 is achieved. Using Eqs. 4 and 5, it is now possible to express f(a/W, D/W) which is sensitive only to variation of t. For doing that only in Eq. 4 all a i s are to be substituted with the functions of Eq. 5. In other words, as t changes, the 24 coefficients b 0 , b 1 ,…, b 23 also changes. Hence, in the next step, for finding the ultimate generalized function, i.e., f(a/W, D/W, t/W), all b i s are plotted against t/W (to get the nondimensional form) and the representative resulting curves are shown in Fig. 7 for the inside case. For outside and average case also, similar curves are obtained. One can be misled by viewing the apparently linear trend of the curves in Fig. 7. This is so as there is very high difference in the order of values of abscissa and ordinate. When regression fitting is made on the data to find dependence of b i s on t/W, it is found that in all cases a second-order polynomial gives a best fit and altogether 72 coefficients (c 0 , c 1 ,…, c 71 ) are derived to correlate b 0 , b 1 ,…, b 23 with t/W with the following representative relations: A regression coefficient of C0.9 is achieved in all cases. Hence, as now all the coefficients c 0 , c 1 ,…, c 71 are known, substituting those values all b i s can be calculated in terms of t/W using Eq. 6 and then using Eq. 5 all a i s can be found in terms of D/W and t/W. And finally using Eq. 4, f(a/W, D/W, t/W) can be obtained. The total general form of the function f(a/W, D/W, t/W) can be expressed as: Equation (7) is general over the ranges of a, D, and t in which the simulation is done and it can predict any shape functional value for all combinations of the above parameters. One should note that FEM gives the freedom of calculating nodal displacement for select nodes specific to requirement and hence it is possible to compare the results across the thickness of the specimen, as has been done in past (Samal et al. 2010a, b) and also in the present case by first individually considering results pertaining to inside and outside nodes and finally with results by taking average of the two cases. In case of experiment, where in past by measurement of compliance, f(a/W) is calculated for select tubes, (Samal et al. 2010a(Samal et al. , b, 2011Sanyal and Samal 2012a, b;Sanyal et al. 2011a, b), it was not possible to find variation of f(a/W) across the specimen wall thickness. Also, through experiment, it is practically not possible to carry out test of specimens covering all the 840 cases as could be done through FEM.

Results
The fitting coefficients for Eq. 7 are shown in Table 2 for all the three cases, i.e., results pertaining to inside and outside nodes and result with average of the two cases. Depending on the need, it is now possible to find expression for f(a/W) for any thin-walled tube, falling within the domain of the tube diameter and wall thickness considered, by proper substitution of values in Eq. 7. The functional form of the fitted functions for all the three cases, i.e., considering results with inside and outside nodes and also with the result taking their average is compared in Figs. 8 and 9 using 3-dimensional plots considering two different sets of combination of the tube geometry. The 2-dimensional projections are also shown within the plots for easy understanding of the comparisons. It is clearly seen that the relative deviation between the results considering the cases inside, outside and average is negligible and for a conservative estimate any of the cases can be considered for calculation of SIF. Also, variation of f(a/W) with change in D/W is much more significant (Fig. 8) than that in t/W (Fig. 9) as can be interpreted from the 2-dimensional projections. In principle, in a general sense it is always better to use the case 'average' for estimation of f(a/W) and SIF to reflect the actual material response as for a real tube through experiment it is impossible to evaluate f(a/W) at select locations across the tube thickness and only a gross estimation is possible. This point is also discussed in detail in Sect. 4. To judge the efficacy of the fitted function, two 3-dimensional plots have been made and shown in Figs. 10 and 11 to compare the output generated through FEM simulation and the fitted function. Considering all individual cases, it is confirmed that in almost all cases, a oneto-one correspondence is achieved and the maximum relative deviation is always \±20 %. For validation of the simulated outcome, the fitted functions are graphically compared with experimental results generated for three different thin-walled tubes with different geometric dimensions used in different reactors (Samal et al. 2010a) in Figs. 12, 13 and 14. The details of the experimental procedure and the analytic derivation for calculation of f(a/W) are available in Ref. (Samal et al. 2010a, b;Sanyal and Samal 2012a) and briefly narrated in Appendix. Considering all the three cases it can be inferred that the relative deviation between the generalized function and experimental result is always \±35 %.

Discussion
From the point of view of classification, there are three broad approaches for calculation of SIF-(1) Analytical, (2) Numerical, and (3) Experimental approaches. Although analytical approaches are the basis for the development of fracture mechanics and they have delivered the basic equations for crack-tip stress and displacement fields, they are not suitable from an engineering point of view as it is not always possible to satisfy all boundary conditions exactly. Thus, in the present study, this procedure is not considered. Since most often, physical processes of interest are lost due to simplifications required for attainment of analytical approximation, numerical methods such as finite pertaining to variation of f(a/W) against a/W and t/W for a fixed D/W of 0.21053 considering results by taking average of simulation results between inside and outside nodes difference (FD), finite element (FE) and boundary element (BE) are extensively used in analyses of fracture problems involving material and geometric non-linearities, for fundamental understanding of the basic physical problem. The availability of high-speed computers and rapid growth in algorithmic methods has helped the wide usage of numerical techniques in fracture analyses. Out of these, FEM allows the analysis of complicated engineering geometries; it enables treatment of three-dimensional problems and it permits the use of elastic-plastic elements to include crack-tip plasticity and thus it is most suitable to handle the current problem numerically. Experimental techniques can never guarantee accurate value of SIF as it is not feasible to get the value directly from experiment.
The variation in size, shape and finish of the notch machined on the specimen for experiment can add incalculable uncertainty in the derived result. In principle, any technique that can measure stresses or displacements can be applied for an experimental determination of the SIF. The most widely applied experimental procedure is the compliance measurement, as used by the authors (Samal et al. 2010a(Samal et al. , b, 2011Sanyal and Samal 2012a, b;Sanyal and Samal 2012a). But the best possible accuracy of compliance measurements can be expected if the points of load application are as close to the crack as possible. In many cases such as the present case, it is not possible to satisfy the above criteria because of the design of specimen and fixture. Hence, although the current FEM results are validated with experimental results, one should note that the large deviation (such as ±35 %) of FEM and fitting results from experimental observation reported does not imply that the prediction made by generalized function is much far from the actual result (Figs. 12,13,14) as experimental result can never be totally accurate unlike analytical methods.
Form Figs. 8, 9, 10 and 11, looking at the 2-dimensional projections, one can infer that there exists approximately an inverse relation between f(a/W) and D/W. In other words, for a fixed t, as D increases, f(a/W) decreases. This fact is interestingly peculiar and investigation for finding a physical basis of this observation may be one of the interesting future scopes of the work. Again, a more interesting observation is noticed if one considers variation of f(a/W) with t/W for a fixed D. For lower D (Fig. 11), there exists a maxima in t/W-f(a/W) plot. Again for moderate D (Fig. 9) and still higher D, a saturation in value of f(a/W) is noticed at higher t/W side. Hence, as D is increasing, the maxima in t/W-f(a/W) plot is shifting towards higher t/W side. Sensitivity of f(a/W) with t/W is also more for higher D. As a future work, this is also very interesting to find out the physical cause of this type of variation of f(a/W) with t. In principle, f(a/W) should be independent of material property and only geometric parameters should determine its value. But in the current context, the same work can be repeated with material having different E value and it can be crosschecked if identical result is achieved. It was also found and highlighted earlier (Samal et al. 2010a) that like f(a/W), g and c functions also depend on the geometry parameters of thinwalled tubes. Hence, a derivation of geometry-dependent generalized g and c functions for calculation of J pl from experimental data is also an interesting future scope of the work.

Conclusions
A generalized shape function expression for calculation of SIF for thin-walled tubular specimens is presented through a detailed FEM modeling and simulation and curve fitting through regression analysis. The function takes care of variation in tube diameter and wall thickness and thus it precludes necessity of carrying out experiment and numerical simulation for tube with any particular geometry within the range and the f(a/W) can thus be predicted with much accuracy. The derived function is validated with available experimental results.
The unusual variation of f(a/W) with D and t is noticed and investigation for finding physical cause behind such observation is suggested with continuation of this work for finding generalized functions for other important functions necessary for calculation of J integral from experiment.
With this method, a more accurate measure of f(a/W) can be found, as the nature of the function is purely analytic and no numerical approximation is associated with the derivation.