Estimation of mixed-mode stress intensity factors using digital image correlation and an interaction integral

. This paper presents a technique for the experimental measurement of stress intensity factors in cracked specimens under mixed-mode loading. This technique is based on full-ﬁeld measurement using digital image correlation and an interaction integral. Such domain-independent integrals are often used in the ﬁnite element method to calculate stress intensity factors. The main advantage of this technique is that the errors made in the estimation of the measured displacement ﬁeld near the crack’s tip do not affect the measurement of the stress intensity factors. The capabilities of the method are illustrated through fracture measurements on compact tension specimens made of marag-ing steel. Another test under mixed-mode loading is presented.


Introduction
Over the past few years, digital image correlation has been developed for the measurement of displacement fields on plane surfaces (Sutton et al., 1983(Sutton et al., , 1986;;Touchal-Mguil, 1997).This technique provides users with full-field measurements and can be used in the determination of fracture characteristics of various materials (McNeill et al., 1987;Anbanto-Bueno and Lambros, 2002) as an alternative to caustics (Rosakis, 1993), Moiré methods (Liu and Ke, 1975) or coherent gradient sensing (Lee et al., 1996).Today, the digital image correlation technique is easy to use and presents the advantage of allowing fast post-processing of the data.At the same time, improvements in the numerical computation of stress intensity factors were achieved in the framework of the finite element method: several techniques were proposed to estimate the energy release rate (Irwin, 1957) or the J -integral (Rice, 1968) and the stress intensity factors by computing domain-independent integrals using the concepts of virtual crack extension (Destuynder et al., 1983;Moran and Shih, 1987) and auxiliary fields (Bui, 1978).These techniques make the determination of local characteristics (stress intensity factors) possible using global quantities (domain-independent integrals).Consequently, such techniques yield accurate and robust results even when significant errors are made in the estimation of the fields close to the crack's tip (Suo and Combescure, 1992).
The main idea of this paper is to combine the two techniques mentioned above: the estimation of stress intensity factors is carried out using domain-independent integrals calculated from experimental data obtained from digital image correlation.This paper focuses on the case of straight cracks under static mixed-mode loading and on the determination of linear elastic fracture mechanics parameters (such as stress intensity factors, fracture toughness or critical energy release rates).Section 1 is dedicated to the use of the digital image correlation technique in ICASOFT (Mguil-Touchal et al., 1997).In Section 2, an interaction integral is developed.Finally, the capabilities of the method are illustrated through several examples.

Digital image correlation
This technique was developed mainly by Sutton et al. (1983Sutton et al. ( , 1986)).Digital image correlation is based on the following principles: the image of the body is described by a discrete function representing the grey level of each pixel.The correlation calculations are carried out for a set of pixels, called a pattern.The displacement field is assumed to be homogeneous inside a pattern.The initial image representing the body before distortion is a discrete function f (x,y) and is transformed into another discrete function f * (x * ,y * ) after distortion or displacement.The theoretical relation between the two discrete functions can be written as where u(x, y) and v(x,y) represent the displacement field for a pattern.Figure 1 introduces some definitions.(For better understanding, the initial and deformed images are represented on the same axis.)The method presented here uses the following homogeneous and bilinear estimated displacement field: u(x, y) The determination of all the terms (rigid body, strains) of this field is based on a cubic spline interpolation of f * and the cross-correlation coefficient: where M is the surface of the pattern of the initial image.Using this technique, the measured displacement field is known at each node of a mesh of the image which consists of pixels with an accuracy of 1/100 pixel.If, in addition, we assume that the displacement field is interpolated linearly between these nodes, this scheme can be considered to be similar to the Finite Element Method (FEM).Hence the idea of using the tools available for the FEM.

Interaction integral
The approach we use is similar to the M-integral proposed by Chen and Shield (1977).Introducing the same concepts, we consider an arbitrary two-dimensional crack.In linear elastic fracture mechanics, the energy release rate in a two-dimensional problem can be obtained from the variation of the Lagrangian of the system associated with a variation of the crack's length.In this section, we present only an outline of the method.Let us consider a cracked domain containing a homogeneous material with linear elastic behavior (λ and µ are Lame constants).The Lagrangian L is defined by the following expression, in which W int is the strain energy and W ext the external work: For mixed-mode separation, we consider a problem with two fields u and u aux .u is the actual displacement field and u aux an auxiliary displacement field.Both fields are kinetically admissible and the Lagrangian of the total displacement field u tot = u + u aux can be decomposed into In the previous equation, L int (u, u aux ) is the interaction Lagrangian.Its expression can be developed as where Different contours for the integration of δL int .
Remarks (i) In the previous equation ,i denotes partial spatial derivatives.
(ii) The external work does not appear in the expression of the interaction Lagrangian because of its linear dependency with respect to the displacement field.
Let us consider a closed two-dimensional domain defined around the crack's tip by its contour ∂ = Ŵ 1 ∪ Ŵ 2 ∪ S + ∪ S − , as shown in Figure 2. Applying the law of conservation of the Lagrangian (Noether, 1918) and assuming that the virtual crack's extension δa is parallel to the tangent at the crack's face on M + or M − and that the crack's faces are traction-free in both the actual and the auxiliary solutions, the following path-independent integral is obtained: where n j are the components of the outward normal vector to the closed domain , and The crack is assumed to be straight.The real field u verifies the momentum balance and the auxiliary field u aux is chosen as the asymptotic crack tip field solution of the problem of an infinite cracked body subject to a tensile stress.Because of the previous assumptions, Q vanishes.If we choose the auxiliary field to be the actual field, we obtain: With these assumptions, is 2J , where J is Rice's J -integral, which is equal to the energy release rate G and can be related to the stress intensity factors by the following Irwin's relation (see Irwin (1957)): In the above equation, K i are the stress intensity factors defined in Irwin (1957).The interaction integral is defined as when the auxiliary fields are asymptotic fields: at the crack's tip, tangent to the crack's faces everywhere. (11) Choosing the virtual crack's extension field q as defined above (see also Figure 3), we can write the interaction integral I int as a domain-independent integral.

Calculation of stress intensity factors
We will now focus on the estimation of stress intensity factors using the interaction integral.In the previous section, the virtual crack's extension field was used to transform the path-independent integral into a domain-independent integral I int .W e also related to Rice's J -integral in the particular case where the auxiliary field are taken to be the real fields.Thus, the J -integral is related to the stress intensity factors through the energy release rate and Irwin's relation (Equation ( 9)).Therefore, here, if K aux i denotes the mode-i stress intensity factor corresponding to the auxiliary field, we can write: Consequently, K 1 (respectively, K 2 ) is calculated by choosing K aux 1 = 1,K aux 2 = 0 (respectively, K aux 1 = 0,K aux 2 = 1).For numerical implementation purposes, the interaction integral is calculated in the domain described in Figure 3. (The area of interest is the domain in which the displacement field is obtained from digital image correlation.)The virtual crack's extension is taken to be x 1 at the crack's tip, and decreases to 0 at the boundary of the domain.For instance: where r is the distance to the crack's tip and R max ,R min are characteristic lengths of the domain.A linear approximation within the elements is chosen based on the experimental data from digital image correlation.Usual finite element numerical integration tools are used to calculate the interaction integral, then the stress intensity factors, from the experimental data.This interaction integral enables one to calculate the stress intensity factors for modes 1 and 2 separately using the measured displacement field from digital image correlation.The procedure is exactly the same for pure mode 1 and for a mixed-mode configuration.Its accuracy is the same in both cases and depends on the accuracy of the measurement of the displacement field.
Remarks (i) Using a virtual crack extension field such as q, the experimental data in the region near the crack's tip (r R min ) are not taken into account in the estimation of the stress intensity factors.Since the interaction integral involves only q derivatives, we avoid zones where the available information on the displacement and strain fields is poor.
(ii) For the estimation of the energy release rate, Irwin's relation is preferred to the J -integral, which can be viewed as a self-interaction integral (the estimated field interacts with itself).Consequently, the error on the measurement of the displacement field perturbs the calculation of the J -integral more than that of the interaction integral, and the error on the estimation of the J -integral estimation is larger.

Mode-1 experiment on a compact tension specimen
The first example concerns a compact tension specimen of the type commonly used for fracture toughness measurements.In this case, the pre-crack is made by electroerosion and, consequently, the radius of the crack's tip (about 0.1 mm) is greater than that of a fatigue pre-crack.The advantage of this specimen's geometry is that an analytical solution was given by Tada for K 1 (see Bui (1978)): where P is the load and a,b,e are the characteristic lengths of the specimen defined in The material is a maraging steel (EZ 2 NKD 18) produced by Aubert et Duval.This material was heated at 480 • C for 4 h.It can be considered brittle compared to other stainless steels, and its nonlinear behavior can be considered perfectly plastic.Consequently, the radius of the plastic zone may extend beyond small-scale yielding assumptions and the measurements of stress intensity factors are likely to be disturbed when the load reaches its critical value (initiation of failure).

Measurement of stress intensity factors
The radius of the area of interest where the displacement field is estimated is 8mm.R min and R max are 3.2 mm and 6.4 mm, respectively (Figure 5). Figure 6 shows the history of the stress intensity factor K 1 for analytical and experimental data.At a low load level, the experimental results and the theoretical solution are in good agreement.When the load reaches its critical value, the values calculated from the experimental data using the method described in the previous section are higher than the theoretical values.This is due to the fact that in this case the small-scale yielding  assumption is not really satisfied (see Section 4.1.2).Finally, K 1 is accurate when the testing conditions follow the assumptions of linear elastic fracture mechanics.
Table 1 shows the sensitivity of the estimation of K 1 to the dimensions R min and R max of the domain.One can see that the results of the calculation of the interaction integral are quite independent of the size of the domain.This illustrates the robustness of the method presented in that the size of the domain has little influence on the results.

Investigation of elastic-plastic behavior
As can be observed on Figure 6, a deviation in the measurement of K 1 occurs for a load value of about 40 kN.It seems that for such loads the small scale yielding assumptions are no longer valid: the size of the plastic zone around the crack's tip increases, and elastic-plastic fields develop and perturb the estimation of elastic stress intensity factors.Let us consider a hardening material with a power-law (Ramberg-Osgood) associated with the uniaxial stress-strain relation:  σ 0 is the reference stress, ε 0 = σ 0 /E the reference strain, α the parameter of the material and n is the hardening exponent.It has been shown that in this case the singularity of the asymptotic field near the crack's tip is modified and Hutchinson-Rice-Rosengren (HRR) type fields develop (Hutchinson, 1968;Rice and Rosengren, 1968).Then, the displacement field is given by In this equation, J is Rice's integral, I n a dimensionless constant which depends on n, and ū a dimensionless function of θ and n.F o rn equal to 1 (respectively, ∞), the material's behavior is linear elastic (respectively, perfectly plastic) and the singularity is r 1/2 (respectively, r 0 ).
To illustrate this purpose, we can study the transition from elastic asymptotic fields to elastic-plastic fields using the experimental value of the displacement field.Figure 7 shows the evolution of the vertical displacement v normalized by √ r.The results, picked along a 45 • radius, are plotted as functions of r/a for different value of P .F o rP = 18 kN, since v/ √ r remains constant in the vicinity of the crack's tip (r/a 0.1), the singularity in √ r is clearly observed.For larger values of r/a, the asymptotic fields do not seem to be dominant.For P = 35 kN, we observe three regions: r/a 0.06, 0.06 r/a 0.125 and r/a 0.125.The second and third regions correspond to those described for P = 18 kN.In the first region, the √ r singularity seems unclear.These three regions evolve as the load increases and for P = 50 kN they correspond approximately to r/a 0.08, 0.08 r/a 0.18 and r/a 0.18.The first region is the largest and the singularity's exponent is much less than 1/2.In the second region, the elastic asymptotic fields are always dominant, but again vanish for large values of r.
For quantitative observations, a plot of the displacement in log-log-scale is shown in Figure 8.In the region of high r/a previously mentioned, the displacement clearly varies like r 1 .Indeed, in this region, asymptotic fields in √ r are superseded by rigidbody rotation displacement and vanish.This observation leads us to an estimate of the region of K-dominance: as P increases, so does this region and at P = 50 kN it tends toward the limit of the area of measurement.At P = 18 kN, elastic asymptotic fields are clearly observed in the vicinity of the crack's tip.At P = 35 kN, the correlation for r/a ≈ 0.06 which can be observed is less good.The transition between these two regions is very clear at P = 50 kN, where the elastic-plastic field develops for r/a ≈ 0.08.This gives us some information on the plastic radius and on the type of singularity which appears when elastic-plastic asymptotic fields develop.At the highest value of P , low-order singularity is observed (near 0, the material is nearly perfectly plastic), which is consistent with the elastic-perfectly plastic behavior of the material.As a consequence, an estimate of the size of the plastic zone at P = 50 kN is r/a ≈ 0.08.Such a plastic radius violates the assumption of small-scale yielding (r p /a 0.02) and explains why the measured values and the analytical values of K 1 differ from each other.
This example also illustrates one of the advantages of using the interaction integral in estimating stress intensity factors.The results shown in Figure 6 were obtained with values of R min and R max corresponding to r/a = 0.1 and r/a = 0.2.At low load levels, where the estimation of the stress intensity factor is accurate, one can observe that the dominance of the K-field in the region 0.1 r/a 0.2 is not wellestablished.However, even in the absence of K-field dominance in the region where the interaction integral is calculated, our technique leads to an accurate estimate of K 1 (or K 1 and K 2 in mixed-mode loading configurations).This is due to the use of asymptotic crack tip fields as auxiliary fields which interact with the experimental field to capture its singularity in the vicinity of the crack's tip.This is the reason why this technique does not require an a priori determination of the region of K-dominance and leads to an accurate measurement of stress intensity factors even when the K-field is predominant only over a very small region around the crack's tip (e.g. at low load levels).
Remark: It has been shown in Rosakis and Ravi-Chandar (1986) that for a given specimen's thickness the three-dimensional zone extends to a radius of about half the thickness.In this work, the three-dimensional zone is a zone around the crack's tip within which two-dimensional plane stress assumption fails.In this zone, the out of plane displacement cannot be derived from the asymptotic stress field using the two-dimensional plane stress constitutive equation.In the results presented here, a telecentric lens was used.Consequently, only in-plane components of the displacement are measured.One can observe, with this technique, that the in-plane components keep a singular behavior within the three-dimensional zone.

Mixed-mode experiment
In order to show that the method is capable of giving the experimental K 1 and K 2 directly, we performed a test in mixed-mode loading configuration using stainless steel tension specimens (120 mm×30 mm×5 mm) with a 40 • pre-crack.The same test was also performed with a 0 • pre-crack.Critical values for pure mode 1 and mixed-mode loading were compared for various fracture parameters: stress intensity factors (K 1 and K 2 ), crack initiation direction (θ c ) and equivalent stress intensity factors (K 1eq and K PM1 ).The first result concerns the crack initiation direction: its postmortem observation (≈−40 • ) is in good agreement with the value predicted from the maximum hoop stress criterion with a critical angle defined by Therefore, the fracture behavior of such a material seems to be governed by the intensity of the tensile stress.An equivalent stress intensity factor K 1eq is given by the maximum hoop stress criterion (see, e.g., Broek (1982)): Another equivalent stress intensity factor, from an energy point of view, is The experimental values of K 1 ,K 2 ,K 1eq ,K PM1 and θ c are presented in Figure 9.As in the previous example, at the highest load level the small-scale yielding assumptions are not satisfied and the behavior of the estimated stress intensity factors is not linear.At low load levels (P 10 kN), the mixing of the modes is not quite stabilized due to measurement errors or placement of the grips.As θ c reaches a stabilized value (P 15 kN), the equivalent stress intensity factors begin to differ.Since the initiation direction is governed by the intensity of the tensile stress, one may assume that K 1eq yields a better estimate of the critical mixed-mode loading level.
Table 2 summarizes the critical values of various fracture parameters for pure Mode 1 and mixed-mode configurations.The critical value of K 1eq for the mixedmode configuration is very close to the pure Mode 1 value.On the contrary, the values of K PM1 are quite different: this energy-equivalent stress intensity factor yields a higher mixed-mode toughness.These results also show that for this material the fracture initiation process is governed by the tensile stress distribution.
Figure 10 shows the distribution of the components of the Green-Lagrange strain tensor around the crack's tip; an arrow points in direction of the predicted crack initiation.The strain tensor was calculated from the displacement field obtained by  digital image correlation.Although digital image correlation seems inaccurate along the crack's faces, these results agree with the fracture mechanics prediction.One can observe that the arrow is an axis of symmetry for these strain distributions, which confirms that the fracture behavior is dominated by tensile stress.

Conclusion
In this paper, first, we briefly presented the digital image correlation technique used for the measurement of displacement and strain fields.Then, we developed a method based on the Lagrangian conservation law for the estimation of stress intensity factors using the concepts of virtual crack extension and auxiliary fields.Neither of these methods is new, but the main idea of this paper, which is to combine both numerical and experimental techniques, is innovative.The first example using CT specimens shows the accuracy and the robustness of the method when the small-scale yielding assumption is satisfied.The measurement of global displacements also provides an interpretation of the transition from the elastic to the elastic-plastic states.
As predicted by elastic-plastic fracture mechanics, a deviation from asymptotic elastic fields to an HRR-type field is observed close to the crack's tip.The second test provides information on the fracture behavior of the material being tested under mixedmode loading.Indeed, the main advantage of the approach presented is that global measurements (displacement fields) yield estimates of local fracture parameters (stress intensity factors) without using perturbed values close to the crack's tip.Thus, it can be viewed as a convenient and promising investigation tool for fracture phenomena.Moreover, since dynamic interaction integrals have already been developed for arbitrarily growing cracks, this method could be used to measure dynamic stress intensity factors successfully.One can also imagine that the asymptotic crack tip fields used as auxiliary fields in the interaction integral could be replaced by high-order singular fields for V-notch specimens, for example, in order to estimate generalized stress intensity factors.The method could also be used for the experimental investigation of various failure criteria, such as the strain energy density.

Figure 1 .
Figure 1.Initial pattern (ABCD centered on P ) and deformed pattern (A * B * C * D * centered on P * ).

Figure 3 .
Figure 3. Domain used to calculate the K i .

Figure 5 .
Figure 5.The area of interest and correlation grid for the CT specimen.

Figure 6 .
Figure 6.Comparison of analytical and experimental values of K 1 .

Figure 8 .
Figure 8. Vertical displacement and power law approximations.

Figure 9 .
Figure 9. Measured values of mixed-mode stress intensity factors.

Figure 10 .
Figure 10.Green-Lagrange strains for the mixed-mode test.

Table 1 .
Influence of the size and position of the domain of integration on the relative error in the estimation of K 1 .

Table 2 .
Critical values of K i (MP a √ m).