Determination of Stress Intensity Factors and J integral based on Digital Image Correlation

Digital image correlation (DIC) is a technique in experimental mechanics to acquire full-field measurement data of displacements and deformations from the surface of specimens or components. Especially for the investigations of cracks it provides additional benefits. The actually present deformation field in the vicinity of the crack tip can be obtained which directly reflects for example crack closure effects or plasticity. Against this background the paper summarizes a procedure to compute the J integral and the stress intensity factors KI and KII based on DIC data. For this purpose the J and interaction integral are computed as line and domain integrals. Through experiments it is shown that the domain integral is less affected by scatter of the DIC data. Furthermore, the specific domain, facet sizes and facet distances slightly influence the results.


INTRODUCTION
racture mechanical experiments are usually conducted and analyzed by the same well-established procedure. Here, a formula for the estimated stress intensity factor for the specimen or structural part is needed. For standard specimens these formulae are found in handbooks [1]. For more complex problems or crack growth in real structures generally simulations are required [2]. Then the experiment is conducted and it is assumed that the actual stress intensity factors (SIFs) are in good agreement with the theoretical solutions. But a lot of external and internal disturbances like small crack deflections, crack closure effects, environmental influences etc. can distort the experiment [3][4][5][6]. One convenient technique to capture the actual loading conditions in terms of displacement and deformation fields is digital image correlation (DIC). DIC is a contactless technique to acquire full field displacements of specimen's surfaces [7]. With the stereoscopic configuration of two cameras 3d information of local displacements can be derived [8]. For this purpose, the DIC software subdivides the recorded images into smaller facets which are arranged like a (self-overlapping) checker pattern. The 3d coordinates of each facet can then be computed based on the pairwise cross correlation of corresponding facets in both stereoscopic photographs. Comparing these 3d coordinates over time finally yields the actual displacement F field on the specimen`s surface during e.g. mechanical loading [9]. The SIFs derived based on this displacement field should consequently reflect the actually acting crack tip loadings. Several authors used DIC to compute fracture mechanical parameters mainly using one of the following two procedures. Firstly, the Williams near field solution can be used to fit its theoretical displacement field (or the corresponding stress field) to the DIC based displacement (or corresponding stress) field (i) [10][11][12][13]. The main advantage of this procedure is that higher order coefficients of the near field and even the out of plane stress intensity factor KIII can directly be obtained. However, one drawback is the basic dependency of the results on the specific fitting procedure. Secondly, the techniques based on integral methods (ii) can be used for the determination of fracture mechanical parameters. Here mainly the J-and the interaction integral J (1,2) are used [14,15]. These integrals can be numerically evaluated as a line integral or an equivalent domain integral (EDI) [16]. One common drawback of all DIC based measurements is the scatter in the obtained displacement fields. As the strain field is the spatial derivative of the displacement field it is particularly sensitive to scatter. But if considering the central limit theorem (CLT) [17] the calculated J-integrals and the corresponding SIFs should be less affected by this scatter as generally multiple discrete data points are used for their calculations. While the interaction integral is just able to determine K I and K II the decomposition method of the J-integral in also able to compute K III based on the DIC data [18]. Generally, digital image correlation is very useful to evaluate the actually present damage process at the (fatigue) crack tip itself [19,20]. Against this background the present paper elucidates important tools and provides a detailed comparison of DIC based fracture mechanical parameters computed using line integrals and EDI.

NUMERICAL PROCEDURES AND IMPLEMENTATION
his section describes how interaction integral and J integral are implemented as a line integral and as an equivalent domain integral (EDI), respectively, into a python-based postprocessor. Both integrals (if properly applied) are path-independent. The interaction integral allows to distinguish into K I and K II . This is not possible based on the J integral. In both cases it is not possible to calculate KIII without any further assumptions. To determine KI and KII the accurate crack tip position is needed. The main advantage of a line integral is its easy implementation compared to the EDI. In return the EDI uses more facet data for a domain of similar size. Especially, this aspect could be a benefit (considering the central limit theorem, CLT) as DIC results are usually more or less noisy. Of course, the same procedures can also be applied to strain and displacement fields computed by finite element simulations. Fig. 1 provides a flowchart illustrating the different steps of the postprocessor. In the following subsection this chart is referenced.

DIC data
As mentioned in the introduction, digital image correlation (DIC) is one method to compute the actual displacement field on the surface of a specimen or component even under complex loading conditions. For 3d DIC all three components of the displacement field are measured, while 2d DIC only captures the in-plane displacement on the specimen's surface. Derived from these data the system is also able to provide the 2d strain tensor on the surface in terms of total strains (i.e. elastic and plastic) (Fig. 1 ①). The computation of the integrals mentioned above also requires the corresponding stress fields. Assuming linear elastic material behavior these can simply be calculated with Hooke's law for plane stress conditions ( Fig. 1 ②) Of course, in this case the integration path or domain must not include plastically deformed regions. Therefore, it is verified facet by facet whether the equivalent stress or strain exceeds the yield criterion (i.e. VM,DIC yield    ). For all further considerations it is assumed that the investigated surface is flat. In general 3d cases the integration procedures must be enhanced to take into account the volumetric distribution of the crack [21].

Line integration
The J integral formulated as line integration is given in Eqn. 2, neglecting body forces and crack-face loads [22]. All of the terms used here are computed from the DIC results.
Here, mn mn is the deformation energy (density), Γ is the integration path, 1 j  = (1 0 0) and j n a normal vector at each path increment ds. Eqn. 3 for the interaction integral   1,2 J is very similar to the J integral, but in this formulation an additional auxiliary field is needed [23]. Here the terms with the superscript (1) are taken from the DIC measurements and the terms with the superscript (2) correspond to the auxiliary field. This auxiliary field is the theoretical near field solution based on the William's series expansion. In this case only the first term of this series is used, as its coefficients are the sought-after values for KI and KII. For that reason the crack tip position must be accurately determined.
Both Eqns. (2) and (3) contain the gradients u i,1 . Here for small displacements u 1,1 is equal to the strain ε xx (definition of coordinate system provided in Fig. 2) and u2,1 needs to be calculated separately. The numerical procedure for the integration is illustrated in Fig. 2. In this illustration the facets are only schematically illustrated and have a more complex distribution in reality. Both integrals are discretized and integrated with the midpoint method [24] [24]. In Fig. 2 (a) the integration path (red line) surrounds the crack tip counterclockwise (i.e. mathematically positive). This path must be outside the plastic zone, as it otherwise loses its path independence. Here the path is aligned parallel or rather perpendicular to the crack. The path is discretized as indicated by the green integration points (Fig. 1 ③). Each integration point contains all values that are required for the integration. For a higher flexibility the position of the integration points does not need to match the facets of the DIC data. The facet data have to be mapped onto the integration points as illustrated in Fig. 2 (b) (Fig. 1 ④). Here the Qhull algorithm [25] is used for triangulation of the nearest three surrounded facets. In the next step linear barycentric interpolation is conducted on each triangle. Both steps are implemented in the python module "scipy.interpolate.griddata" (SciPy 1.1.0, 2018). For the calculations of u 2,1 (= Δv / Δx) additional auxiliary points are needed as also shown in Fig. 2 (b). These additional points are generated by shifting the integration points in the positive and negative x direction. Then for each facet the difference quotient is calculated using Eqn. 4 ( Fig. 1 ⑤).

Equivalent domain integral
The following section gives a detailed description of the methods to evaluate the J and interaction integrals as a domain integral. In this case the formulae for the J integral (Eq. 5) and interaction integral (Eq. 6) are almost similar to Eqns. 2 and 3 [26]. But instead of the normal vector n j, a q-function is needed which will be described in more detail later. Furthermore, the domain formulation contains an integration over an area increment dA instead of a line increment ds.
The next step is the discretization of area A for the subsequent numerical integration. The finite element method provides appropriate tools. Fig. 3 illustrates this procedure. The integration domain surrounding the crack tip in subfigure (a) is covered with a mesh of elements. In this formulation linear QUAD4 elements with four nodes are used, as they provide sufficient accuracy ( Fig. 1 ⑧). The green lines around the mesh indicate the value of the q-function. At the inner border it has to take the value 1.0 and at the outer border 0.0 [15]. Similar to the line integral, the crack tip and the crack path itself must be excluded. The following aspects must be taken into account when defining the domain. Firstly, as a kind of inherent artefact of the DIC procedure the crack faces themselves generally show unrealistically high strain values, as the crack path itself is not automatically excluded from the calculations. Secondly, (and as mentioned above) it must be taken into account that Hooke's law is not applicable inside the plastic zone. Fig. 3 (b) shows exemplarily one local QUAD4 element. Again, the nodes do not have to match with the positions of the facets. The mapping of the facet data onto the nodes is conducted with the same linear interpolation algorithm as described above for the line integral ( Fig. 1 ⑨).
For the numerical integration with the Gaussian quadrature the next step is the mapping of the nodal data (blue) to the integration points (green) (Fig. 1 ⑩). This second mapping is conducted with the (linear) ansatz functions of the QUAD4 element [27]. The ansatz functions are formulated in the local coordinate system (ξ, η) of the element (compare Fig. 3 (b)). Then, any value can be interpolated from the nodes to any position inside the element. In this case these values are coordinates, displacements, stresses and strains. For the J and interaction integral also the gradients of the displacements ui,1 are needed. The gradients are calculated by using the Jacobian determinant. For the Gaussian quadrature a four-point integration scheme is chosen. The q-function can be interpreted as a function describing a virtual crack extension in the crack propagation direction [16]. As mentioned above, it takes values between 1.0 on the inner border and 0.0 at the outer border as illustrated in Fig.  3 (a). To fulfill these requirements the domain is subdivided into several smaller subdomains as shown in Fig. 4 (a) (Fig. 1  ⑪). Each subdomain requires an individual formula as summarized in Eqn. 7. This is necessary as there is no closed formula for this type of domain. The q-function is plotted in Fig. 4 (b). It has constant transitions between the different subdomains. But this contour leads to discontinuity in the corners of C L and C R in the derivative functions shown in Fig.  4 (c) and (d). Due to the rectangular shape of the domain these discontinuities are unavoidable. This shape was preferred instead of a circular shape as it is more flexible and much more space of the image section could be used. Furthermore, the q-function should not significantly influence the results of the J or interaction integral.
Finally, all values required for solving the J and interaction integral are available ( Fig. 1 ⑬). As mentioned before, for the computation of the interaction integral an auxiliary field is needed ( Fig. 1 ⑥, ⑫), as given in the next part (Eq. 8).

Auxiliary field
As mentioned above, the computation of the interaction integral requires an auxiliary field. All values with the superscript (2) correspond to this field; please see Eqns. (3) and (6). For the determination of K I and K II the first term of the Williams series expansion is utilized. Generally, the Williams field describes the stress field in the vicinity of the crack tip under linear elastic conditions. Here, K I and K II represent the coefficients of the first term. The corresponding stresses and displacements are summarized in Eqn. 8 [26] [26]. The elastic strains can be calculated from the stresses using Hooke's law (see Eqn. 1). I  I I  11  11  11  2  2  2  I  I I  I  I I  22  22  22   I  I I  2  12  12  12 2 2 I  I I  2  2  1  1  1  I  I I  I  I I  2  2  2  2   2 2 2 2

Calculation of stress intensity factors
The final step is the calculation of the stress intensity factors (Fig. 1 ⑭). Two calculations (2a) and (2b) have to be conducted with the interaction integral using different coefficients for the auxiliary field. In the simplified Eqn. 9 the first calculation is executed with    are used [26].
The restriction in this system of equations is that alternately KI or KII of the auxiliary field are zero. In the numerical implementation both calculations can be simply conducted in the same loop. The numerical code used in this work has been verified using results of finite element simulations. Except small deviations due to numerical inaccuracies the line and domain integral lead to identical results for pure mode I, II and their mixed mode scenarios. Furthermore, under linear-elastic conditions the J integral can be transferred into terms of K with Eqn. 10 which can be interpreted as some kind of representative or equivalent stress intensity factor under mixed-mode conditions. Under pure mode I loadings KJ is equal to KI.

EXPERIMENTAL RESULTS AND DISCUSSION
he following section presents selected results for both the line and domain integration based on experimental DIC data.

Experimental setup
The mixed-mode device in Fig. 5 is a special loading frame (a) designed for testing of mixed-mode (single or combined mode I and II) conditions, i.e. different KI/KII ratios [28]. For a rotating angle of β = 0° as shown in the picture the specimen is subjected to pure mode I loadings, while rotating this device to β = 90° leads to pure mode II loadings. A sketch of the so called compact-tension-shear (=CTS) specimen is shown in (b). For further details including the calculations of the stress intensity factors please refer to [19] and [28]. In this work, the specimen has overall dimensions of 90 mm x 145 mm and a thickness of 5 mm. Because of its welldocumented mechanical properties, the aluminum alloy AA2024-T3 was used which has a Young's modulus E of 71 GPa and a yield strength σy of 350 MPa. The mixed-mode device is attached to a standard servo hydraulic testing machine with a load capacity of 60 kN. Like in the picture the crack started from the left side using a saw cut as crack initiation site (acut  40 mm). The specimen was cyclically loaded with ΔK = 10 MPa√m and R = 0.  ). The reflex camera has the advantages that the setup is comparable cheap but still providing sufficiently high optical resolution. In this setup with one camera only 2d (in plane) measurements are possible. Consequently, the distance between camera and specimen must be kept constant during testing, and the specimen surface must remain plane. For a high spatial resolution the lens was operated at the shortest possible distance. The corresponding image section was 21.9 × 14.6 mm. Compared to the 12M setup a spatial resolution of 0.06 mm or 250 Pixel/mm has been achieved with a facet size of 19 x 19 Pixel and a facet distance of 15 Pixel.

Experimental results
For the sake of brevity only selected experimental results are presented in the following section, which provide an insight into the main characteristics of both integration procedures. In all cases the exact crack tip position was determined by measuring the crack manually after the experiment to avoid inaccuracies in the computation of KI and KII. Fig. 6 and Fig. 7 show results of the line and domain integration technique for pure mode I loadings and mixed-mode loadings. Usually, a pure KI scenario has the highest significance in technical applications for straight crack propagation, as mixed-mode loadings generally lead to kinking cracks. While Fig. 6 was obtained by the GOM stereo camera system from the rear side, Fig. 7 was obtained by the reflex camera at the front side. Both analyses have different spatial resolutions as described above. The subpictures (b) and (c) show the local results of the domain integration in terms of KElement (i.e. element-wise) and illustrate how the scatter of the DIC measurements affects this domain. An equivalent graphic for the line integral is shown in subpictures (e) and (f) for the corresponding integration points. In this case each line in subfigure (e) and (f) represents a different integration path whose results are shown as data points in subfigure (d). For domain and line integration it is important that no node or integration point is inside the plastic zone because the stresses are calculated with Hooke's law which is only valid in the linear-elastic case. Furthermore, the DIC system cannot accurately compute strains close to the crack path. Therefore, the strains and derived stresses in the vicinity of the crack are unrealistically high which is exemplarily shown in the subfigure (a). Consequently, these regions must intentionally be excluded from the integration procedure. For the pure mode I load case in Fig. 6 the subfigures (b, c) and (e, f) show that the color distribution of line and domain integral are very similar at first glance. It has to be mentioned that the plots should not be compared quantitatively due to the different number of integration points used. Qualitatively the plots (b) and (d) reveal a symmetric distribution against the y-axis for K * Element. The highest values (red) are located at x = 0 mm which is directly above and below the crack tip. Here the equivalent stress as shown in subfigure (a) has the highest amount. shows the distribution of the results for K I and K II for different integration paths (red curves). The individual values (red squares) oscillate around their mean value (green line) which basically indicates that their variation is mainly due to the inevitable scatter in the DIC results. Systematic deviations would occur if the DIC results were erroneous, the crack tip position would be incorrect or a part of the integration path would be located within the crack path or plastic zone. The spread for KI in this example is about 1.3 MPa√m which is only ~4.5 % from the corresponding mean value. Similar results are shown for KII with a spread of 0.9 MPa√m. These observations clearly indicate that in case of line integrals multiple integration paths should be computed as it is highly likely that only one path would not reliably represent the local crack tip loadings. For the mode I case in Fig. 6 values for KJ based on Eqn. 10 are KJ,area = 27.53 MPa√m and KJ,line = 28.28 MPa√m. Here the difference between both integration procedures is significantly lower compared to KI calculated with the interaction integral. Therefore, this value can also be used for comparison of the results. Fig. 7 shows a mixed mode case where the load frame was rotated by 45 deg.  Both analyses of the pure mode I and the mixed-mode case revealed that for the mode I and II different regions of the integration domain yield very different amounts of local K * I,Element and K * II,Element values. While for K I the regions above and below the crack tip have the highest contributions, in case of K II all areas around the crack tip are crucial because of the locally negative and positive values. Therefore, KII should be basically less sensitive to local (stochastic) scatter as larger regions of the DIC-field significantly contribute to its computation. The local scatter that can be seen for directly adjoining integrations points is more an effect of the comparable small values. Furthermore it can be seen that the distribution for KI is more smoothly.   In the next step the influence of the facet sizes and distances are investigated based on the domain integral. As shown above, this integral is less vulnerable to general scatter and should therefore reveal in how far the facet size affects the computation of KI. Fig. 9 summarizes the results intentionally using very large facets for the sake of clarity. These large facets with a size of up to 90 Pixel and a distance of up to 45 Pixel (see subfigure (a)) basically enhance the accuracy of the measurements as in general the precision of the calculated displacements is higher. In return, the local resolution is lower compared to the smaller facets as evident in subfigure (b) showing much more details of the characteristics of the crack tip field. With a facet size of 19 Pixel and a facet distance of 15 Pixel even small details like the shape of the plastic zone at the crack tip can be revealed. But comparing (a) and (b) the scatter in the computed stress field in (b) is clearly higher than that in (a). For that reason, it was exemplarily investigated how facet size and distance influence the results of the integration procedures. Fig. 9 (c) shows the results for the domain integral for facet sizes between 19 and 90 Pixel and facet distances between 15 and 45 Pixel. The bars show that the results are almost independent of the facets and the local scatter. The deviation of the highest and lowest value is only < 1 %. Furthermore, the element size of the mesh was also varied and set to 0.5 mm (648 Elements) and 1.5 mm (72 Elements). But again, number and size of used elements do not significantly influence the results in this example case. Therefore, it is concluded, that the results are mainly dominated by the overall quality of the picture acquisition for the DIC procedure and only to a lower degree by the details of the subsequent analyses. This conclusion should be kept in mind especially if only one camera is used not allowing to capture the potential out-of-plane deformations during the experiments. This would result in an erroneous displacement field which in turn would result in computing K-factors that could significantly deviate from the actually prevailing local loading conditions. Finally, the influence of the crack tip position itself was investigated (see Fig. 10) as it is needed for the calculations of the auxiliary field in the case of the interaction integral. The example case shown in Fig. 6 was used for the variation of the crack tip position under mode I loading. To avoid disturbances caused by variation of the integration paths because of scatter of the DIC data all seven paths were kept constant. During the computations the crack tip position was virtually shifted within a region of 2 x 2 mm 2 around its actual location (here x = 0.0, y = 0.0) as indicated by the colored squares in Fig. 10. This should also reflect the achievable range of accuracy if the real crack tip position is estimated manually based on the contour plot of the plastic zone. The results for K I in Fig. 10 (a) reveal that the deviation between the minimum and maximum values is just about 0.18 MPa√m. This error is less than 1 % compared with the nominal value. Along the crack direction (x-axis) the values computed for KI are almost constant while they increase when shifting the crack tip in the positive y-direction. For KII the range is just 0.36 MPa√m which is also negligible if compared to the deviations discussed in the previous chapters. Here K II decreases in the crack direction while it remains almost constant for shifts perpendicular to the crack. Summarized, the results for the mode I case show that the crack tip position itself does not significantly influence the results compared to the scatter of the DIC data (see Fig. 6 (d)).

SUMMARY AND CONCLUSIONS
n algorithm for the analysis of mode I and II stress intensity factors and J integral based on interaction and line integral is presented in detail in the first part of the paper. These integrals were implemented as a line and an equivalent domain integral. The main purpose is to derive KI, KII and J directly from experimental results based on digital image correlation analyses. Therefore, experiments with a mixed-mode load frame were conducted. The main conclusions can be summarized as follows: (1) The implementation of the J and interaction integral as a line integral is easier compared to the domain integral, because no separate element formulation is needed. (2) Because of the same type of the J and interaction integral both integrals can be computed in the same loop in the line and domain integral program code. (3) The combination of DIC and the J and interaction integral provides reliable data for the experimental determination of K I and K II . (4) For the line integral several integration paths must be analyzed to obtain reliable data. (5) The domain integral is clearly less affected by scatter of the DIC results as generally many more integration points are used compared to the line integral. (6) The computation of the J integral does not require the exact spatial coordinates of the crack tip; consequently, the J integral basically provides more robust results than the interaction integral.