Extraction of stress intensity factors for 3D small fatigue cracks using digital volume correlation and X-ray tomography

This paper describes a methodology used to compute Stress Intensity Factor values along the curved front of a fatigue crack inside a nodular cast iron. An artificial defect is introduced at the surface of a small sample. The initiation and growth of a fatigue crack from this defect during constant amplitude cycling is monitored in situ by laboratory x-ray tomography. The method for processing the 3D images in order to compute SIF values is described in detail. The results obtained show variations of the stress intensity factor values along the crack front.


Introduction
In the last decade, in the transport industry, a large research effort has been devoted to the reduction of the structural weight of vehicles in order to reduce fuel consumption and comply with tougher environmental regulation. Globally, this has lead to a more severe use of materials which are submitted to higher mechanical loads than in the past. For theoretical models which aim at predicting the durability of materials this is a challenging evolution. In the field of cyclic mechanical behaviour (fatigue), the Paris law established in the sixties is widely used to model fatigue life of components based on propagation laws.
Such laws are established with standard test procedures which monitor crack propagation rates of a relatively long (generally a few millimetres) through crack which front is assumed to be straight. In practice, however, it might take a very substantial part of the fatigue life of a component for it to contain a fatigue crack comparable in size and shape. While the size issue can be addressed via the similarity principle (up to some extent, see for example (author?) [1] for a discussion on that topic), the effect of the 3D shape of the crack and the variation of stress intensity factor values that it implies is generally neglected.
To establish crack propagation laws of fatigue cracks with a curved front in opaque materials, several techniques have been proposed which are either direct (beach marking) or indirect (compliance measurements) and which all suffer from various drawbacks (see ref [2] for a discussion). The availability of powerful x-ray beams which can go through several millimetres of metals, has made it possible to analyse directly the mechanics of crack tips in the interior of an opaque material either by direct imaging or by analysing diffraction patterns or by combining both signals [3].
The emergence of digital correlation methods, on the other hand, gives access to crack tip displacement fields from which it is possible to extract Stress Intensity Factor (SIF) values. This has been done at the surface [4,5,6,7] and also in the bulk [8,9,10] of cracked samples. In the latter case, the images used are obtained by x-ray tomography and the need for a high spatial resolution restricts the size of the studied samples to a few millimetres. In practice, the samples studied in 3D so far all contain crack tips extracted from larger coupons. The stress state at the tip of those cracks is therefore largely unknown because it is inherited from a larger crack and also because it is likely to have been affected by the machining of the small sample.
This paper describes a methodology used to overcome these issues. An artificial defect (notch) is introduced at the surface of a miniature metallic fatigue sample. The initiation and growth of a fatigue crack from this defect during constant amplitude cycling is monitored in situ by laboratory x-ray tomography. The processing of the 3D images required to compute SIF values from the recorded 3D images is described in detail. The potential of the method for computing SIF values along the front of curved crack fronts is demonstrated.

Material properties
The material used in this study is a nodular graphite cast iron (3.4 wt.% C, 2.6 wt.% Si, 0.05 wt.% Mg, 0.19 wt.% Mn, 0.005 wt.% S, and 0.01 wt.% P) with a homogeneous dispersion of graphite nodules which are used as natural markers for Digital Volume Correlation (DVC) [11]. Casting and appropriate heat treatments produce a ferritic matrix microstructure (grain size ∼ 50µm) with a 14% volume fraction of graphite nodules (average diameter 45 µm). The static and cyclic mechanical properties of the material are listed in Ta  The aim of this work is to characterise cracks initiated and grown in situ. In this material, within smooth samples, fatigue cracks usually initiate on pores which are produced by the solidification process. Such defects are relatively rare. For various reasons [12], the size of the fatigue samples characterised by tomography samples has to be kept small (e.g. section 1.6 × 1.6 mm 2 , gage length ∼ 5mm in this work). The probability to find such a casting defect in a sample extracted randomly from a bar is, hence, very low. An alternative solution has been adopted: cracks are initiated from well controlled artificial sharp notches. In a previous work, FIB machining had been used to produce such notches in an Al sample [2]. For the material studied here, however, the notches that could be produced with this technique in a reasonable time (a few hours) were too shallow (maximum depth ∼ 20µm) to act as an efficient initiation site. The notches were then produced by laser machining [13] using a 1 kHz repetition rate regeneratively amplified Ti:Sapphire laser system employing chirped pulse amplification. A pulse length (full-width half-maximum) of approximately 40 fs is measured using an auto-correlator. A femtosecond laser is chosen as it produces a negligible heat affected zone around the notch due to the high intensity of the laser pulse which directly vaporizes the material. For all experiments presented here, rotating polarization obtained with a spinning zero order half wave plate was used. The beam is focused with a 10x microscope objective using an energy per pulse of 50 µJ. Cast iron samples are placed on a computer controlled x-y translation stage for precise sample positioning under the laser beam. The notches are produced by translating the stage back and forth under the laser beam until the desired depth is obtained.
The results shown in this paper have been obtained with a notch located at the centre of one of the lateral surface of the sample (surface crack) 2 . The notch is perpendicular to the applied cyclic stress and has approximately a square shape 180 × 180 µm 2 (see figure 1) with an initial opening of 15µm and a tip radius (both measured by optical inspection at the sample surface) of about 7µm.

X-ray tomography
A laboratory tomograph has been used. The X-ray source is operated at an accelerating voltage of 85 kV with a transmission tungsten target. The minimum source/detector distance imposed by the machine (see below) is 15 mm resulting in a voxel size of 2.5 µm. A current of 174 µA and a mean of 3 radiographs with an exposure time of 1 second have been used; 720 angular positions have been taken which corresponded to a scan duration of 45 minutes. The reconstruction time (weighted filtered back projection, Feldkamp algorithm) on a Desktop computer with a bi Intel xeon quad core processor, 3 GHz CPU, 64 GBytes RAM, takes approximately 15 minutes for a 560 × 560 × 900 voxels 3D image.

In situ cycling
Cyclic loading of the sample is performed in situ using a miniature hydraulic fatigue machine [14]. In order to decrease the voxel size in the reconstructed 3D images, the sample must be as close as possible to the laboratory X-ray source (divergent beam configuration), thus the sample is placed into a long slender aluminium tube which diameter and thickness have to be kept as small as possible to be transparent to X-rays. In this work the diameter of the tube (Al 2024 alloy) is 16 mm and the thickness 1 mm [15].
As the contrast of the crack on the radiographs is blurred by the presence of numerous graphite nodules, crack growth in the interior of the sample is assessed by recording a tomographic scan. The crack shape evolution is studied through the residual of the DVC which correspond to the difference between the reference and the deformed and corrected volumes as described in [11].
Once the sample is set up in the loading stage, two scans are recorded at low (40 MPa) and high stress (290 MPa) to characterise the initial state. The sample is then cycled in situ at a constant stress amplitude of 300 MPa (R=0.1). One hundred thousand cycles are applied prior to the first recorded scan. Then 2 scans are recorded every 15 kcycles, one at 40 MPa and another at 290 MPa, to monitor the crack evolution as described in [16]. Because fatigue cycling is performed at constant stress amplitude, the crack accelerates during propagation; near the end of the experiment, smaller steps are used between two scans according to the approximate crack speed da dN observed on the residuals.

Optical flow equation and resolution
The optical flow equation assumes that the grey level variations are due to the displacement of the material points only. Its generic form is written as f (X) = g(X + u(X)).
In this equation, X is the vector denoting the position of the current point with respect to the image or volume coordinate system and u the unknown displacement vector field. f and g, are the reference and deformed grey level functions representing the image or volume of the sample. Once a basis functions has been chosen to describe the variations of the displacement field in space a non-linear least-squares resolution strategy is adopted. The unknow vector U is then obtained as a minimizer of The minimization of η 2 results in the resolution of a series of linear systems until convergence is reached. More details on the resolution strategy can be found in the references listed in the next sections.

Displacement measurements
For the estimation of the displacement field, a finite element basis function is adopted as proposed in (author?) [17] for the 2D case and in (author?) [18] for the 3D case. This basis function is enriched in the spirit of the eXtended Finite Element Method [19] by using discontinuous enrichment functions. Discontinuous enrichment functions of the Heaviside type are adopted to measure the displacement discontinuity across the crack faces as in (author?) [20] which is dedicated to 2D measurements. In the present paper, this strategy is extended in 3D.
Due to the noise in the data acquisition chain (X-ray source, scintillator, CCD) and to the ill-posed nature of the problem (a vector field is searched for whereas only a scalar equation can be written), resolution/uncertainty compromise must be solved. Previous studies have shown that the uncertainty which alter the measured displacements increases when the element size decreases [17,18,20,21]. Thus, when trying to capture crack tip field, which exhibits strong variations i.e. , requires high spatial resolution, small element size is needed. Without any further improvement, the required element size would lead to very high noise level. To circumvent this drawback, a regularization is added to the DVC equation. For the analysis of a fatigue cycle, the small scale yielding assumption is fulfilled and an elastic behaviour of the material can be adopted within a loading cycle. Based on the work by (author?) [22] in 2D, an elastic regularization is added to the minimization problem. In the formalism then introduced by (author?) [23,24], a cut-off wave length c below which the variations of the displacement field are governed by linear elasticity, is introduced. This regularization uses the equilibrium gap method and it enforces that the linear elastic balance of momentum equations are satisfied. Using this strategy, and contrary for example to particle tracking techniques, there is no need to have a sufficiently high marker density. The weakness of the image texture is compensated for by the mechanical regularization that enforces smoothness and elastic balance of momentum constraint. The regularization includes the traction free boundary condition on the crack faces and also on the sample surface. This kind of boundary conditions holds meaningful mechanical information and it is thus important for the finite element mesh to incorporate both an accurate description of the crack faces and of the sample surface. For these purposes, the DVC formalism is also used to extract an accurate geometrical description of the sample boundaries and of the crack surface.

Contour and crack detection
From the acquired grey level volumes, the optical flow equation is used to compute the 3D displacement field between two scans. As the sub-surface depth of the crack analysed in this work remains small (∼ 500µm), displacement fields must be computed very close to the sample surface, specially in the first stages of growth. Because an elastic reularization is used, the description of stress free surfaces must also be accurate. This is achieved by using the optical flow equation as explained in the following paragraphs. A similar strategy is also used for determining the crack surface.
The detection of these surface relies on the same problem which is to extract the silhouette of an object in a 3D image. Due to the graphite nodules, a direct thresholding of the x-ray scans is neither adapted for extracting the crack surface nor for meshing the sample surface. For the crack surface, the correlation residual (a voxel reconstruction of the optical flow residual) can be used as proposed in (author?) [20] or (author?) [25]. However due to the small voxel size used, direct thresholding of this residual volume is again not accurate enough. A methodology based on the optical flow equation is used as follows. Inspired by (author?) [26] and (author?) [27], the conservation of grey levels is written between a virtual scan f * and the reference one f .
f * is a virtual volume that describes a silhouette of simple geometrical shape with an ideal grey level transition between black and white. b is a grey level function used to account for the lower spectrum of the actual scan. s is the displacement of the virtual object surface so that the actual reference scans f at the modified surface location matches the virtual scan. For the sample surface detection, the virtual volume shape is a cylinder and for the crack surface, a virtual plane is used. The variations of s are described using B-splines functions and the same non-linear least-squares minimization as for the displacement measurement is used.

Stress intensity factors calculation
In previous 2D analyses [28,29,30,31], Williams' series [32] were extensively used to extract stress intensity factors and crack tip position directly from digital images [33] or by least-squares projection [34] of the displacement field. This strategy has also been used in 3D for through cracks with nearly straight crack fronts by (author?) [9,10]. In the present experiment where curvilinear crack fronts are studied, the original methodology needs some improvements.

Williams' series
Considering an infinite homogeneous elastic solid with isotropic behaviour, Williams [32] has obtained the in-plane and out-of-plane solution for the stress distribution produced by a semi-infinite crack. The solution is decomposed over a series of displacement functions,indexed by n for each mode m, u n m and v n m , along the crack direction and orthogonal to the crack direction respectively. These in-plane fields for mode I and mode II are written in a polar coordinate system (r, θ) as u n I (r, θ) = r n/2 κ cos and u n II (r, θ) = −r n/2 κ sin v n II (r, θ) = r n/2 κ cos For mode III, only the out-of-plane displacement w n III is described: w n III (r, θ) = r n/2 cos(n(θ − π)/2).
Based on a linearization of the perturbation due to a mispositioning of the crack tip, it was shown by (author?) [33] that the crack tip position can be estimated by using the first mode I super-singular field (n = −1). An iterative search is performed so that the position of the equivalent elastic crack tip, i.e. corresponding to a vanishing amplitude of the super-singular term, is progressively reached. This iterative search is initialized with a semi-elliptical crack front extracted from the DVC residuals as shown in Figure 4d. Then at each iteration of the procedure, the displacement interpolated on circular projection domains is projected onto a basis containing displacement functions for n varying from −3 to 9 for the 3 modes. As long as the equivalent elastic crack tip position is not reached, for a given position along the crack front, the first mode I super singular amplitude is not vanishing and the stress intensity factor estimation given by n = 1 amplitudes is not reliable. Once a sufficiently small shift (typically less than a tenth of voxel) to the equivalent elastic crack tip is estimated, then the procedure is stopped and the stress intensity factors values and crack tip positions are stored for this position along the initial estimate of the crack front.

Adaptation to 3D crack front
First, the projection onto the Williams series is performed over circular domains which orientation depends on the local normal to the crack surface and to the crack front. As illustrated on Figure 2, the projection domain is adapted to the local crack geometry in order to reduce the approximation due to a misalignement of the local crack tip coordinate system with the projection plane. Note that during the iterative search of the crack tip position for a given position along the front, the projection domain is moved along the local normal to the initial semi-elliptical crack front. The second improvement consists in using the modification suggested by (author?) [35] to account for crack curvature. It is shown in [35] that the expression of the stress field in the crack front vicinity is modified when considering the perturbation due to locally curved crack morphology. From these corrected expressions of the stress field, a corrected displacement basis is obtained. This corrected basis is then used to extract the stress intensity factors and crack tip position from the measured displacement field by a least-square projection. The method has been validated by tracking curved crack fronts within displacement fields obtained from linear elastic finite element simulations. Robust and accurate detection is obtained when such noise free displacement fields are used.

Results and analysis
In this section, a shallow surface crack of maximum depth around 500 µm is analysed for three different numbers of fatigue cycles (120, 140 and 160 kcycles). Figure 3 shows a displacement field obtained using the proposed methodology. Note that in order to avoid the elastic regularization prescribing the crack front position, a small cylindrical region over the crack is front is excluded from this regularization for which the cut-off wave length was set to 256 voxels. The mesh size in the vicinity of the crack front is 8 voxels.
As mentioned before, elastic regularization of the DVC problem requires an accurate description of the crack and sample surfaces; Figure 4 shows the results obtained using the method described in section 3.3 for the crack surface after 160 000 fatigue cycles. It can be seen on this figure that the method used gives a good description of the slightly curved sample surface and of the crack surface topography. About the latter, the method used to extract the crack surface tends to smooth the crack deflections which are usually observed in this type of material when graphite nodules are encountered by the crack tip [36]. However as shown on Fig.4c, some roughness of the crack surface with a typical range of ±12 voxel i.e. ±30 µm is detected; this value is close to the graphite nodules average size.
The detection of the crack front using the method described in section 4.1 is illustrated in figure 5 for the crack observed after 160 000 fatigue cycles. The corresponding distribution of K I values along the crack front is shown in figure 6 for three different crack positions corresponding to different numbers of fatigue cycles. These results have been obtained by projecting the displacement field measured between two scans acquired at the minimum and maximum loads of the same loading cycle. The extraction is performed using discs of 100 voxels in diameter and Williams functions for n from −3 to 9. The region in the vicinity of the crack front which is not submitted to elastic regularization is not used for the extraction as well as a band of 40 voxels along the crack faces.
To the author's knowledge, this is the first time that K values are extracted by from experimentally measured displacements along the front of semi elliptical low fatigue cracks initiated and propagating in situ from a notch in an optically opaque material. Compared to previous results published in the literature, in the same material [37,11] or in Al alloys [8,38] studied by particle tracking, the cracks studied here have built their own plastic zone. This one has not been altered by the preparation of the sample from a larger specimen. If the displacement fields studied in this work are therefore more realistic, the orders of magnitude of the crack openings are however at the limit of what can be achieved by the correlation method (see for example the much larger values of crack opening values in [16] for a comparison) which has to be "assisted" by the introduction of an elastic regularization. At the intersection between the crack front and the specimen surface, the results given by the method might be considered with some care due to lower robustness of the crack front detection. This can be due to the use of the Williams' series that may not be appropriate for describing the mechanical state in the vicinity of the crack vertex.
As mentioned before there is no equivalent experimental data available to compare with our results. X-ray diffraction which has been used by several authors for studying crack tip displacements [39,40] in the bulk of metals for through cracks does not have the necessary spatial resolution for the curved crack fronts of the short cracks studied here. One can however compare the experimental SIF curves of figure 6 with calculations either analytical or making use of the finite element (FE) method. We have shown elsewhere [41] that the Raju and Newman analytical approximation [42] gives SIF values which are very close to FE calculations (less than 5% discrepancy) even in the case of a curved crack front showing some irregularities (elastic calculation), therefore in this paper our experimental results are only compared with the analytical calculation ( figure 6).
This comparison is qualitatively satisfying: the maximum (resp. minimum) SIF values are measured at the surface (resp. in the bulk) of the sample (140 and 160 kcycles curves) in agreement with the trend observed on the analytical curves. The latter being a good average of the measured ones for the three crack fronts investigated. When looking at the variation of the SIF values along the crack fronts, it can be seen that the range of variation between the maximum (surface) and minimum (bulk) values is much larger for the experimental curves (∼ 9 MPa √ m) than for the analytical ones (∼ 1 MPa √ m). The size of the cracks detected by the Williams series at the surface (2c) is smaller than the size of the crack detected in the bulk (a), this gives the characteristics U shape of the K I curve shown on figure 6. In other words, the crack front is lagging behind at the surface. This so-called tunnelling effect which is well known for fatigue cracks in macroscopic CT specimens cycled with a load ratio close to 0.1 is usually attributed to a larger closure effect at the surface (plane stress -larger plastic zone) [1]. It has also been observed for elliptical cracks [1] and corner cracks [43] and even in the case of small (millimetric) tomographic samples [2].
The values of K shown on figure 6 are close to the typical value of K op which can be found in the literature for the cast iron studied (between 4 and 6 MPa √ m [11]) and therefore such closure effects are likely to happen. The fact that, in nodular cast iron, closure plays a role during the propagation of such very small cracks had already been suggested by Clément et al. [44]. Our results tend to show, in addition, that the level of closure might vary along the front of the crack.
Finally, figure 7 shows the comparison with the experimental crack growth rates measured in this work with some literature data on the same material [45]. The points obtained at the surface and in the bulk in this work fall in the characteristic short crack zone of the Paris curve; the size of the crack at 160 000 cycles being probably at the transition between the short and long crack regime. This behaviour is again similar to that observed by by Clément et al. [44] for artificial cracks of comparable sizes. Besides, the points measured in this work at the surface tend to fall closer to the nominal curve while bulk measurements are closer to the effective curve. This is in agreement with the above mentioned argument of a lower level of closure in the bulk.
Compared to previously published results [11,16], the experimental setup has been improved in order to avoid the effect of an initial unknown plastic state. One drawback, however, is that x-ray tomography limits the size of the sample. Therefore the physical size of the crack is obviously limited. As a result, graphite nodules that are used for supporting the DVC analysis introduce strong heterogeneities in the local material properties resulting in strong interactions between crack propagation and the material microstructure. This can be seen for example on figure 5 where the irregularities of the crack front detected by the Williams series can be correlated with the presence of a porosity and /or nodules. As discussed above, the SIF values which were extracted from the measured displacement fields are nevertheless coherent with simulations and results on large samples. It can be expected, however, that a material with smaller microstructural markers (but large enough to be imaged by the tomography) would increase the quality of the results obtained.

Conclusions
Laboratory X-ray tomography has been used to monitor the initiation and growth of a fatigue crack in a nodular cast iron sample. SIF values along the curved crack front were extracted from displacement fields by DVC. Compared to previous work where the tips of through cracks coming from larger coupons were examined, this method has the advantage of producing a crack with its "own" plastic zone/stress field . A laser machined notch is used to initiate in situ the small fatigue crack (physical size of the order of 500 µm). Some tools have been developed in 3D to map with a reasonable accuracy the sample and crack surfaces automatically. The small size of the studied crack (< 500µm) and the low values of the crack opening which had to be measured are at the limit of of what can be detected with classical correlation methods for describing crack tip displacement fields with a high enough spatial resolution. This problem was solved by assisting correlation with a mechanical regularisation. The SIF values obtained are in agreement with analytical calculations and reveal a typical short crack behaviour. Larger values of SIF were found close to the sample surface; this was attributed to a higher level of closure compared to the bulk of the specimen. The method presented here is limited by the relatively large size of the nodules, which act as microstructural markers for the correlation, compared to the physically small size of the crack.

Acknowledgements
The measurements described here were carried out under ANR funding granted through ANR-09-BLAN-0009-01 (RUPXCUBE Project). The authors would like to acknowledge Dr. M.C.Baïetto, Prof. A.Gravouil, Dr. F.Hild, Mr. C.Roux and Dr. S.Roux for fruitful discussions during the PhD work of Joël Lachambre. Figure 3: a)Example of a displacement field for a surface crack after 160 000 fatigue cycles (Mode I opening -stress direction along z). The shape of the crack in the bulk of the sample is also shown. The "smooth" part of the field corresponds to the part of the sample which is regularized whereas displacements within the "speckled" regions close to the crack tips at the sample surface (arrows) are only driven by DVC. b) Map of the crack opening.  Figure 5: Superposition of the detected crack front position and the material microstructure. The colormap shows the estimated distance to the crack front using the first supersingular term of the Williams' series (see the text for details). The black line indicates the iso-0 line of this distance, i.e. the detected crack front. The laser notch from which the crack is initiated is also shown as well as graphite nodules and a pore.