Controlling Stress Intensity Factors During a Fatigue Crack Propagation Using Digital Image Correlation and a Load Shedding Procedure

Digital Image Correlation is used for controlling load shedding fatigue crack propagation. A specific algorithm is used to perform Stress Intensity Factors (SIFs) and crack length estimation in real time. Crack length measurements are validated by comparison with potential drop technique. SIFs results are compared with more common techniques using standard analytical formula considering confined plasticity at the crack tip. The proposed non-contact method is shown to be a powerful tool to control crack propagation.


Introduction
Investigating fatigue crack propagation laws requires Stress Intensity Factors (SIFs) and crack length (a) estimation as functions of the number of cycles (N) [1,2]. Those parameters need to be estimated accurately during the fatigue test for a robust identification of the fatigue law parameters. A large number of measurement techniques exists to estimate crack length during its propagation. Each method can be chosen depending on its accuracy and sensor constrains. Compliance techniques can be used to perform crack length measurements in real time [3]. Potential drop [4] is also an accurate method which has already been used to detect short crack during fatigue propagation [5]. On the other hand Crack Opening Displacement or Crack Tip Opening Angle are carried out to study crack propagation [6,7]. From a given crack length estimation, SIFs are then calculated using standard analytical formula depending on the material elastic properties and the sample geometry [8]. Such procedures need a very accurate calibration that is often specified for normalized specimens or that requires a calibration procedure for specific specimen geometries. Moreover those measurement tools are based on global data (loading force) and they can be sensitive to boundary conditions applied on specimens. Therefore, whereas their robustness was demonstrated in some standard cases, they can be poorly adapted when experimental conditions are constraining. For example, performing SIFs measurements into aggressive or high temperature media becomes complicated as the use of measurement sensors in direct contact with the specimen is often compromised. In these cases, it is more convenient to use non-contact techniques [9]. For Stress Corrosion Cracking studies e.g., potential drop technique could disturb electrochemical measurements [10].
As other alternative non-contact methods such as photoelasticity [11] or thermoelasticity [12], Digital Image Correlation (DIC) provides the surface in-plane displacement field of a loaded sample. It is an accurate tool to identify fracture parameters such as both crack length and SIFs [13][14][15][16]. During fatigue crack propagation, plasticity appears at crack tip, which can induce local crack closure. These non-linearities can also affect SIFs measurements. Irwin suggested a method to correct the SIFs estimation based on crack length reestimation taking into account plasticity at crack tip [17,18]. An alternative route consists in using a specific basis functions to represent the displacement field around the crack tip such as Williams' series with non-singular and super-singular terms. The latter allows considering plasticity influence on SIFs measurements [19].
On the other hand, fatigue crack propagation laws are often identified using the load shedding method which consists in modifying the loading path during the propagation [20,21]. This method can be used to estimate K I threshold or to impose a decrease of the Fracture Process Zone (FPZ) size [22]. Contrary to a propagation with a constant loading path [23], stability of the crack propagation can be obtained with this technique. It enables a more accurate study of the crack velocity as a function of the stress intensity factor. Moreover an uncontrolled increase of K I involves growth of plasticity in the vicinity of the crack tip which may affect the accuracy of fracture parameters estimation.
The present work is aimed at developing a real time SIFs and crack length estimation procedure using DIC in order to control the crack propagation conditions. First the load shedding procedure is presented. Then, the fracture parameters identification using real time DIC is exposed. The performances of the proposed technique are then illustrated by several sensitivity analysis. Lastly, controlled fatigue crack propagation results using the load shedding method are presented and compared with results obtained with the usual combination of potential drop technique and the use of analytical formulas.

Fatigue Test and Load Shedding Procedure
A MTS hydraulic tensile machine is used to apply a loading path with an initial "fatigue" sinusoidal profile F 0 (t) (see Fig. 1): (1) The hydraulic system allows reaching quite high frequencies, up to f = 20 Hz. The loading ratio R is set to 0.1 as recommended by standards defining the production of pre-cracks [23]. The maximal value for the force during one cycle for the initial profile is F 0 max . Load shedding can be used to predict the stress intensity factor threshold [2,21]. It is also used to control the crack propagation to obtain for instance a constant size of the fracture process zone r p which is estimated in plane stress by: where K Imax is the maximum value of the mode I stress intensity factor during the cycle and σ y the initial yield stress of the material. Mode I stress intensity factor K I is usually estimated as : A geometrical correcting function F g (a, w) depending on the crack length a and the sample width w allows accounting for finite size effect [8]. F max e w √ πa gives the value for a semi-infinite crack in an infinite medium. In this formula, e is the thickness of the sample and F max the maximum force applied during the actual cycle. With a constant loading path, the increase of a induces therefore the concomitant increase of the stress intensity factor. The consequence is both the increase of r p and a misleading estimation of K Imax using the analytical formula when r p becomes too large. In order to control the size of the fracture process zone, it appears therefore important to control the value of K Imax .
A real time estimation of K Imax is then required to use a load shedding procedure. In the following, digital images of the sample surface are acquired and a DICbased method for a real time extraction of K Imax is proposed. Image acquisition and DIC are performed every 100 cycles during the test to control K Imax . We consider, that between two image acquisitions, the crack increment is small. Then, we assume the relative variation of F g (a, w) √ πa as negligible compared to those of F max and K Imax . Thus, a simple relation is obtained between the relative variations of F max and K Imax : this relation being valid under the assumption of small increments of a. Using this formula a step-by-step regulation of F max can be developed based on the estimation of dK Imax and K Imax without using any analytical formula as equation (3). We adopt a step-by-step (a step being a group of 100 cycles) regulation to impose the maximum value of the mode I stress intensity factor K n+1 Imax at step n to a target value K t Imax . At stage n, we have By recurrence, the loading at step n + 1 is computed using the following equation: where the maximum load to apply (F n+1 max ) is given as a function of both the initial maximum load F 0 max and the ratios between the target K t Imax value and estimated values of K k Imax for the previous steps. The regulation is applied using a feedback procedure. An initial loading path is applied on the sample (F 0 (t) in equation (1)). Then every 100 cycles, image acquisition and DIC are proceeded to extract K n Imax . Thus, at each step n, F n (t) is deduced from α n calculated from all previous image acquisitions and equation (1): The real time DIC procedure is described in the next section. It enables to estimate α n which is converted onto an analogical signal to be used as a feedback parameter to impose the loading path. Figure 2 illustrates the architecture of the load shedding experimental device with DIC.

Displacement Field Decomposition
Williams' series enable to describe the displacement field around a semi-infinite straight crack into a linear isotropic elastic material [24]. The displacement field is expressed in the complex plane as an infinite series: where j is the loading mode (I or I I), and n is the order of the corresponding term in the series : n = 0 accounts for rigid body translations, n = 2 corresponds to rigid body rotation and T-stress, n = 1 gives directly the standard Westergaard asymptotic displacement field around the crack tip which is directly proportional to SIFs.

Amplitudes Estimation from Digital Images
Digital images are used in order to extract the amplitudes corresponding to the above described decomposition of the displacement field. DIC is a full-field displacement measurement technique based on the optical flow equation that states for the conservation of the grey level between two different pictures of the surface of a loaded specimen. The displacement field u is searched for between a reference digital (grey level) image, f and a deformed one, g. A non-linear least-squares technique is adopted. It consists in finding u that minimizes being the region of interest. A finite number of n j is chosen for n between n min and n max . If l is a generic indice l = (n − n min + 1) + ( j − 1)(n max − n min + 1), the solution, i.e. the set of amplitudes ω l , is obtained iteratively by solving a sequence of linear system where and To reduce computational cost, ∇g x + u is usually replaced by ∇ f x in the expressions above. This allows computing M kl and also part of F k once for all. As ∇g x + u = ∇ f x when convergence is reached, this modification only affects the convergence rate of the algorithm. The amplitudes of the Williams' series functions are thus directly extracted from digital images as proposed in [15]. Note that these functions depend on the distance to the crack tip and also on the angle with respect to the crack orientation. Therefore, crack tip location and crack orientation must be estimated accurately. Considering mode I fatigue crack propagation, the crack orientation is easy to obtain with a really low level of uncertainty. Conversely, estimating the position of the crack tip along the crack line is more difficult. Based on [15,16], the following recurrence formula is derived: Considering that the crack tip was mispositioned by dx along the real axis (see Fig. 3), the decomposition of u on the basis functions calculated using the shifted coordinates is given by whereω n j are the amplitudes of the shifted n j functions. Note that for the shifted basis n < 0 terms are also considered. Indeed, using equation (15), a 1st order Taylor development of equation (16) and identifying the corresponding terms, the following relation is obtained: crack tip θ r x y dx origin of basis functions  For n = −1, the mispositioning of the crack tip is thus written as equation (18) In References [15,25], the position of the crack tip is estimated through an iterative procedure (Fig. 4): For the next fatigue cycles, one also has to estimate the load shedding factor in order to maintain a constant stress intensity factor. The latter has to be estimated at each cycle even for a non-zero mispositioning of the crack tip. By considering n = 1 in equation (17), a linear relation betweenω 1 j ,ω 1 j ,ω 3 j and dx is obtained: K I is then estimated by K I = 2μ √ 2πω 1 j where μ is the shear modulus of the material. Thanks to the linearization of equation (16), one can estimate the actual crack tip position and stress intensity factors even if the basis functions used for DIC do not correspond to the actual tip position. Of course the linearization is only valid within a limited range of dx. Depending on this interval, a number of basis functions must be generated a priori. The next section is devoted to the estimation of the range of confidence and also to the performances of the proposed real time DIC extraction procedure for the crack length and mode I SIF. The number of order and the size of the Region Of Interest have to be optimized to improve the computational cost. The influence of these parameters on the performances of the method as well as the validity of the linearization are investigated in the next section.

Fatigue Test Set Up
Fatigue tests are performed on Zircaloy (Zry-4) samples which are notched to localize crack initiation. The sample is designed with a section S = 2.5 × 8 mm 2 . The notch is 0.5 mm in length. During the fatigue test a 1,200 × 1,600 pixels digital camera is used to take pictures of the sample surface. A speckle is applied on it to generate a random grey level distribution (see Fig. 6). It is applied using a white then a black paint spray. With this technique the stain size order of magnitude is about ten pixels. A telecentric objective is use to avoid artificial dilatation due to off-plane displacement. To avoid transient effects the image acquisition is driven without stopping the tensile machine with a very short exposure time (dT = 500 μs on Fig. 1). A high powerful lighting is applied allowing the use of such a small exposure time. To warranty good reproductibility, pictures are taken for a given value of the loading force (F meas ). It was chosen to acquire images at a constant ratio F meas /F max equal to 0.8 (see Fig. 1). The image acquisition is thus synchronized with the loading force. The image acquisition and digital image correlation are performed while the test in running and fracture parameters are estimated using the procedure presented in the previous section.

Validity of the Linearization
Calculations with different sets of basis functions ( n j (x, y)) shifted of dx up to ±20 pixels from the known crack tip position were made on a picture with a defined crack and K I magnitude. Thanks to the method proposed above, ω n j are deduced fromω n j . The estimated crack shift normalized with the width of the Region Of Interest (w r on Fig. 6), dx/w r is given on Fig. 7(a). A very good agreement is obtained between the applied shift and the one estimated by DIC. The relative error in K I estimation (ε K I ) is shown in Fig. 7(b). It is shown that this error increases with the applied shift. However, for dx/w r = ±0.03, the relative error remains lower than 0.03%. The ratio dx/w r = ±0.02, gives a relative error lower than 0.1%. Thus at least 1 0.02 = 50 basis functions have to be used to ensure an interpolation error on K I extraction lower than 0.1%.

Influence of the Domain Size on K I
The influence of the Region of Interest (ROI) size is also investigated to find an optimum between accuracy and computational time. For this study, the number of terms taken in the Williams' series expansion is fixed at 7 and −3 for respectively positive and negative indices of n. The horizontal position and the width of the ROI (w r ) is fixed with the width of the specimen (650 pixels in our case) (Fig. 6). The vertical position of the ROI is chosen to center the notch on it. Several calculations are performed by increasing the ROI height (h r ) for different deformed pictures with different known crack lengths (a/w) and K I magnitudes. Figure 8 gives the evolution of the normalized quantity, ε hr = K I (h r )/min| h r (K I ) which is the ratio between K I magnitude for a given value of h r and the minimum value of K I for a range of h r (min| h r (K I )). For each picture, Fig. 8 . 8 Influence of the zone of interest size on the K I estimation that the number of pixels to get stabilized values for ε hr increases with the crack length. Starting from low ROI size, ε hr decreases down to its converging value as more and more information is available to constrain the solution. After a critical ROI size, ε hr deviates from its stabilized value. One may invoke that more terms in the Williams' series should be incorporated in order to account for far-field effects.
To combine accuracy and fast calculation, we suggest to impose ε hr ≤ 1.02 which ensure an error lower than 2%. In this configuration, Fig. 8 reveals that with w r = 650, h r has to be larger than approximatively 600 pixels.

Influence of the Number of Terms on K I
Super-singular terms (n < 0) and non-singular terms (n > 1) are required to perform a correct SIFs estimation. Non-singular terms are dominant far from the crack tip and they allow capturing far-field and finite size effects. The maximum number of terms is discussed herein to find the minimum number of terms in order to ensure robustness and accuracy. For a deformed picture with a given crack length (a/w = 0.5), K I extractions are performed using different number of non-singular terms for different number of super-singular terms. The latter is defined by n min , the number non-singular terms being n max . The study is performed with −3 ≤ n min ≤ −1 and 1 ≤ n max ≤ 12. Figure 9 shows that whatever the  Fig. 9 Influence of the Williams' series order number on the estimation of K I number of super-singular terms n min , at least seven nonsingular terms are needed to obtain a stabilized value of K I . Further it is shown that n min = −3 is required in order to get stabilized influence of super-singular terms. Crack tip non-linearities can be described by super-singular terms as it has been demonstrated by Henninger et al. [19]. Consequently they need to be considered into the Williams' series expansion to give an accurate description of the displacement field around the crack tip.
The same study is performed onto different crack lengths for the same magnitude of K I (1/3 ≤ a/w ≤ 2/3). ε K I defines the relative error in K n I extraction, considering n max non-singular modes into the calculation compared to the calculation with n max = 12(K 12 I ): Figure 10 shows the evolution of ε K I as a function of n max and evidences the importance of considering high order terms when the crack is longer. Indeed ε K I remains almost negligible when a/w 1/3 for n max ≥ 5, although it gives more significant error magnitudes when the crack becomes longer (∼ 10%) for n max = 5. Thus at least seven non-singular terms have to be High order non-singular terms allow correcting finite size effects. These effects are more important when the crack length increases because the distance between the crack tip and the stress free edge of the specimen decreases. Thus for long cracks, the displacement field is more sensible to boundary effects and more high order terms have to be added in the displacement field decomposition.
It can be illustrated by the displacement maps along y direction generated by different n I functions as shown in Fig. 11. On the one hand, super-singular terms are more sensitive in the crack tip vicinity (Fig. 11(a) shows the −3 I field). On the other hand, high order of non-singular terms govern far field and finite size (close to the boundaries) effects (Fig. 11(c) shows the 7 I field). Figure 11(b) gives the displacement map generated by 1 I which gives the K I estimation. Last, the resulting displacement field considering all modes between −3 and 12 is depicted on Fig. 11(d).

Crack Propagation Control: Results
The presented method is applied on the Zry4 sample with a chosen target value of K t Imax . To avoid undesirable effect of crack initiation we only focus on the crack propagation. A pre-crack is initiated until a/w = 0.25. The load shedding procedure is then applied by controlling K I as described before with real time estimation of K I using Williams' series. Contrary to constant loading propagation, constant K I propagation ensures a constant crack propagation velocity. The regulation procedure enables to make the crack propagate until a/w = 3/4. Figure 13 shows how the K I estimation evolves during the propagation and how the target SIF (K t Imax ) value is accurately maintained. Yet, DIC gives a K I estimation using kinematic field that can differ from the analytical formula based on elasticity when confined plasticity develops. For comparison purposes, potential drop technique is used to estimate the crack length and a semi-analytical formula is derived for K I estimation.

Crack Length Evolution
Potential drop (Direct Current Potential drop-DCPD) is a conventional method to estimate crack length during propagation. Crack length is estimated by measuring the potential drop on each side of the crack while the specimen is run through by a constant electric current. The crack length is deduced using Johnson's formula [4]. Figure 12 shows the comparison between the crack length obtained by real time DIC and potential drop. A very good agreement is obtained that validates crack length estimation with the proposed real time DIC.

Stress Intensity Factor Control
Usually, SIFs estimations are obtained using analytical formula which are derived from elasticity theory by introducing a geometrical correcting function (equa-tion (3)). For an edge crack, the geometrical function is given as The formula is given for an applied stress. In our experimental device, a specific vice grip has been Moreover the material behaviour plays an important role. Plasticity at the crack tip could affect SIFs estimation. Irwin et al. [17,18] suggest a method to take those effects into consideration. Corrected SIF (K * Imax ) is deduced by considering that the crack length is incremented by r p obtained from equation (2). The corrected SIF K * Imax is thus given by K * Imax =F * g (a, W) · σ π a + r p . Figure 13 shows the comparison between the DIC using Williams' series K I extraction (solid line) and the classical analytical estimation in dotted line (K * Imax ) (the crack length being estimated with potential drop technique). Similar results are obtained. This procedure enables to ensure a constant fracture process zone size (r p = 0.13 mm). This confirms the performance of the proposed technique as well as its potential for running experiments under controlled SIFs.

Conclusion and Perspectives
A specific DIC technique algorithm using the Williams' series formalism has been developed for a real time estimation of K I and crack tip position. The method was successively compared with potential drop technique. It was used to assist the crack propagation by controlling the loading intensity during the propagation. It enables to get more stabilized crack propagation and therefore to perform more accurate SIFs estimation along a more widespread field of propagation. The proposed method allows investigating crack propagation laws for a given K I magnitude and crack growth rate and can be used to define easily Paris propagation laws. The non-contact properties of the method allow its use in aggressive media where measurements based on sensor are commonly complicated to run. Interesting applications of this method are expected in the field of Stress Corrosion Cracking studies.