An Image Processing Technique to Identify Crack Tip Position and Automate Fracture Parameter Extraction Using DIC: Application to Dynamic Fracture

Abstract

Background

A challenge for experimental fracture mechanics studies using vision-based methods is the accuracy with which the crack tip can be located in the region of interest for extracting fracture parameters. When using full-field displacement measurement methods such as digital image correlation (DIC), positioning the crack tip coordinate system could greatly influence the accuracy of stress intensity factors for brittle materials.

Objective

The objective of the present work is to develop improved methods of tracking crack tip position for fracture parameter extraction for problems involving moving fracture fronts (e.g. dynamic crack growth).

Methods

An improved image processing-based automated method for identifying the location of a propagating crack tip is proposed here. The primary inputs to the method are two-dimensional displacement fields measured using DIC. An edge detection methodology using a series of partial derivative computations is used to locate the crack tip.

Results

The proposed method’s performance is verified using simulated displacement fields with a sequence of controlled crack tip positions for mode I and mixed-mode examples. The method is used to locate crack tip positions from mixed-mode dynamic fracture experiments and extract instantaneous stress intensity factor histories. Consistency is shown between baseline and automated methods and post-initiation stress intensity factor histories varied by approximately 5% with the maximum variation being under 10% for the mixed-mode experiments.

Conclusions

The automated fracture parameter extraction method produced consistent results with those extracted using traditionally accepted methods, indicating that the proposed automated approach is a marked improvement due to its systematic nature and processing efficiency.

Appendices

Appendix A

Non-maximum Suppression Procedure

A simple example with an edge running in the horizontal or vertical direction can be used to illustrate the concept of non-maximum suppression. A 5 × 5 excerpt from a random magnitude intensity field with an edge oriented in the vertical direction is shown in Fig. 18 with resulting simple gradient computations.

Fig. 18

Example intensity gradient a with edge running in vertical direction, and b example resulting gradient computation

The gradient magnitude is calculated from the individual directional gradient values using equation (5) and the directions of the gradient vector can be computed using equation (6).

With the magnitude and direction of the partial derivatives known, the edge points can next be separated from the non-edge points. This information is used to adjudicate points within the field of gradient values that may be an actual edge, with the objective being to arrive at an edge that is exactly 1 data point wide. The first step here is to organize the partial derivative values according to direction, such that they are grouped into bins, [0°, 45°, 90°, 135°]. For the example problem illustrated here, the angles are tabulated in Fig. 19. For this example, since the values are primarily dominated by the vertically oriented edge, all of the directions round to 0°.

Fig. 19

Gradient direction calculations (left) and direction values binned to the nearest 45° increment (right)

With the directions known, each value can then be compared to the eight data points that surround it. More specifically, each point is compared to its neighboring points only in the direction of the angle of the gradient. For instance, if the direction is determined to be closest to the 45° direction, the data point is compared to the point to its upper right and lower left. The value at the given data point is then taken as the maximum of the 3 points along that direction. In the current example, the direction values are all 0°, thus each gradient value is only compared to its left or right neighbor as illustrated in Fig. 20.

Fig. 20

Original gradient values (left) with arrows showing the general direction along which maximum values are determined, resulting in the thinned matrix (right)

Appendix B

Stress Intensity Factor Extraction Using Over-deterministic Least-squares Approach

For the displacement field around the crack tip prior to crack initiation, the over-deterministic least-squares results can be computed using the equations reported in [39] for the crack sliding (u_x) and crack opening displacements (u_y):

$$\begin{array}{l}{u}_{x}=\sum\limits_{n=1}^{N}\frac{{\left({K}_{I}\right)}_{n}}{2\mu }\frac{{r}^\frac{n}{2}}{\sqrt{2\pi }}\left\{\upkappa\cos\frac{n}{2}\theta -\frac{n}{2}\cos\left(\frac{n}{2}-2\right)\theta \right.\\\qquad\quad \left.+\left\{\frac{n}{2}-{\left(-1\right)}^{n}\right\}\cos\frac{n}{2}\theta \right\}+\sum\limits_{n=1}^{N}\frac{{\left({K}_{II}\right)}_{n}}{2\mu }\frac{{r}^\frac{n}{2}}{\sqrt{2\pi }}\\\qquad\quad \left\{\upkappa\sin\frac{n}{2}\theta -\frac{n}{2}\sin\left(\frac{n}{2}-2\right)\theta +\left\{\frac{n}{2}-{\left(-1\right)}^{n}\right\}\sin\frac{n}{2}\theta \right\}\end{array}$$

(12)

$$\begin{array}{l}{u}_{y}=\sum\limits_{n=1}^{N}\frac{{\left({K}_{I}\right)}_{n}}{2\mu }\frac{{r}^\frac{n}{2}}{\sqrt{2\pi }}\left\{\upkappa\sin\frac{n}{2}\theta +\frac{n}{2}\sin\left(\frac{n}{2}-2\right)\theta \right.\\\qquad\quad \left.-\left\{\frac{n}{2}+{\left(-1\right)}^{n}\right\}\sin\frac{n}{2}\theta \right\}+\sum\limits_{n=1}^{N}\frac{{\left({K}_{II}\right)}_{n}}{2\mu }\frac{{r}^\frac{n}{2}}{\sqrt{2\pi}}\\\qquad\quad \left\{-\upkappa\cos\frac{n}{2}\theta -\frac{n}{2}\cos\left(\frac{n}{2}-2\right)\theta +\left\{\frac{n}{2}-{\left(-1\right)}^{n}\right\}\cos\frac{n}{2}\theta \right\}\end{array}$$

(13)

In the preceding equations, µ is the material shear modulus, and r and θ are the polar coordinates with crack tip as the origin and $\kappa =\frac{3-\upsilon }{1+\upsilon }$ for plane stress. The coefficients K_I and K_II, when n = 1, are the mode I and mode II stress intensity factors. For digital image correlation experiments, the u_x and u_y fields are known for a set of points in the polar coordinates r and θ.

By selecting a group of points in the vicinity of the crack, a set of equations can be formed to determine coefficients ${\left({K}_{I}\right)}_{n}$ and ${\left({K}_{II}\right)}_{n}$. Using an over-deterministic approach, the experimental crack opening displacement can be used for extracting mode I fracture components whereas the crack sliding displacements can be used for mode II fracture components. However, it has been shown that by transforming experimentally measured in-plane Cartesian displacements into radial (u_r) and angular (u_θ) components, more accurate SIFs can be found in mixed mode problems [7]. That is, the Cartesian displacement components can be transformed into polar components as shown in equation (14).

$$\left\{\begin{array}{c}{u}_{r}\\ {u}_{\theta }\end{array}\right\}=\left[\begin{array}{cc}\cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{array}\right]\left\{\begin{array}{c}{u}_{x}\\ {u}_{y}\end{array}\right\}$$

(14)

For the analyses presented in the present work, the radial (u_r) components are utilized for computing the SIFs. Using this technique, these equations can be expanded out to any number of higher order terms. For the present work, the equations were expanded for up to 10 terms and stress intensity factors were taken once K_I and K_II were determined to have converged. Measured displacement data was extracted for $0.5\le {^r}/{_B} \le 1.5$ and $-120^\circ \le \theta \le 120^\circ$. The over-determined equation set was formed and solved for minimizing the least-squares error to compute values of K_I, K_II for the crack up to the point of initiation at a range of values of n.

Once the crack begins to propagate, the opening and sliding displacements can instead be written as:

$$\begin{array}{l}{u}_{x}=\sum\limits_{n=1}^{N}\frac{{\left({K}_{I}\right)}_{n}{B}_{I}\left(C\right)}{2\mu }\sqrt{\frac{2}{\pi }}\left(n+1\right)\\\qquad \left\{{r}_{1}^{n/2}\cos\frac{n}{2}{\theta }_{1}-h\left(n\right){r}_{2}^{n/2}\cos\frac{n}{2}{\theta }_{2}\right\}\\\\\qquad +\sum\limits_{n=1}^{N}\frac{{\left({K}_{II}\right)}_{n}{B}_{II}\left(C\right)}{2\mu }\sqrt{\frac{2}{\pi }}\left(n+1\right)\\\\\qquad \left\{{r}_{1}^{n/2}\cos\frac{n}{2}{\theta }_{1}-h\left(\overline{n }\right){r}_{2}^{n/2}\cos\frac{n}{2}{\theta }_{2}\right\}\end{array}$$

(15)

$$\begin{array}{l}{u}_{y}=\sum\limits_{n=1}^{N}\frac{{\left({K}_{I}\right)}_{n}{B}_{I}\left(C\right)}{2\mu }\sqrt{\frac{2}{\pi }}\left(n+1\right)\\\qquad \left\{{-{\beta }_{1}r}_{1}^{n/2}\sin\frac{n}{2}{\theta }_{1}-\frac{h\left(n\right)}{{\beta }_{2}}{r}_{2}^{n/2}\sin\frac{n}{2}{\theta }_{2}\right\}\\\qquad + \sum\limits_{n=1}^{N}\frac{{\left({K}_{II}\right)}_{n}{B}_{II}\left(C\right)}{2\mu }\sqrt{\frac{2}{\pi }}\left(n+1\right)\\\qquad \left\{{{\beta }_{1}r}_{1}^{n/2}\cos\frac{n}{2}{\theta }_{1}+\frac{h\left(\overline{n}\right)}{{\beta }_{2}}{r}_{2}^{n/2}\cos\frac{n}{2}{\theta }_{2}\right\}\end{array}$$

(16)

In the above equations, µ is the material shear modulus, and r and θ are the polar coordinates with crack tip as the origin as before and $\kappa =\frac{3-\upsilon }{1+\upsilon }$ for plane stress. The longitudinal and shear wave speeds are defined as ${C}_{L}=\sqrt{\frac{\left(\kappa +1\right)\mu }{\left(\kappa -1\right)\rho }}$ and ${C}_{S}=\sqrt{\frac{\mu }{\rho }}$, respectively. The non-dimensional quantities, ${\beta }_{1}=\sqrt{1-{\left(\frac{c}{{C}_{L}}\right)}^{2}}$ and ${\beta }_{2}=\sqrt{1-{\left(\frac{c}{{C}_{S}}\right)}^{2}}$ are used to compute the spatial variations of ${r}_{m}=\sqrt{{X}^{2}+{\beta }_{m}^{2}{Y}^{2}}$ and ${\theta }_{m}={\mathrm{tan}}^{-1}\left(\frac{{\beta }_{m}Y}{X}\right)$ based on the crack speed, c. Also, B_I, B_II, D, and h are defined in equation (17).

$$\begin{array}{c}\begin{array}{cc}{B}_{I}\left(c\right)=\frac{\left(1+{\beta }_{2}^{2}\right)}{D}, & {B}_{II}\left(c\right)=\frac{2{\beta }_{2}}{D}\end{array}\\ D=4{\beta }_{1}{\beta }_{2}-{\left(1+{\beta }_{2}^{2}\right)}^{2}\\ h\left(n\right)=\left\{\begin{array}{cc}\frac{2{\beta }_{1}{\beta }_{2}}{1+{\beta }_{2}^{2}}& \mathrm{for\ even}\ n\\ \frac{1+{\beta }_{2}^{2}}{2}& \mathrm{for\ odd}\ n\end{array}\right.\\ h\left(\overline{n }\right)=h\left(n+1\right)\end{array}$$

(17)

Appendix C

Mapping from FE-space to Uniformly Gridded Space

As previously stated, the FE model is comprised of a variety of element shapes due to geometry around the crack, and therefore, does not create output on a uniformly spaced grid of points. A mapping procedure is thus necessary to create a uniformly spaced grid of displacement data for testing the algorithm. To that end, an inverse FE mapping technique was created to generate displacement fields. In the mapper developed for the current effort, the input file for the source finite element model contains all of the node numbers, nodal coordinates, and element connectivity. For each of the 4-noded elements in the finite element model, the mapper locates the grid points that reside within its boundaries using a polygon search algorithm coded in MATLAB^®. Since the element shape could be in the form of any four-sided polygon, potentially distorted, a numerical routine was then used to determine the parametric coordinates of each of the destination grid points within the space of their parent source element in the original model. An example of the dissimilarity between the two data point locations is shown in Fig. 21. A set of x- and y-coordinates on a uniform grid was created at the desired “output” point locations, as shown in red. The source model elements and nodes are shown in black. The relationship between the global space and the parametric space is also illustrated in Fig. 21 with an example map-to point shown by the dark-shaded point, x_p.

Fig. 21

Example to show dissimilarity between simulated data and gridded data near crack faces (top) and mapping from global coordinate space to parametric coordinate space (bottom)

The global coordinate of any point within the boundary of the element is a function of the parametric equation, N, and the global coordinates of the nodes that define the boundary of the polygon. For a 4-noded quadrilateral element, the global coordinate of a point, x_p and y_p,, is defined in equation (18).

$$({x}_{p};{y}_{p})=\sum_{i=1}^{4}{N}_{i}\left(\xi ,\eta \right)({x}_{i};{y}_{i})$$

(18)

where i is the node number, x_i and y_i are the global coordinates of the i-th node and ξ and η are the parametric coordinates.

The parametric equations, N, for a quadrilateral element are [40]:

$$\begin{array}{c}{N}_{1}=\frac{1}{4}\left(1-\xi \right)\left(1-\eta \right)\\ {N}_{2}=\frac{1}{4}\left(1+\xi \right)\left(1-\eta \right)\\ {N}_{3}=\frac{1}{4}\left(1+\xi \right)\left(1+\eta \right)\\ {N}_{4}=\frac{1}{4}\left(1-\xi \right)\left(1+\eta \right)\end{array}$$

(19)

The parametric coordinates, ξ and η, for the target output point can be located using an iterative procedure. For a given iteration, the parametric space is split up into a 5 × 5 grid of points. The values of ξ and η, are used to calculate the resulting global coordinates at each of these points on the 5 × 5 grid. The point within the grid that results in coordinates that have the shortest Euclidean distance to the actual point of interest is used as the initial guess of the next iteration. That initial guess becomes the center point of a smaller 5 × 5 grid that is part of a subdivision of the grid in the previous iteration. This iterative process continues to subdivide the parametric space into smaller and smaller 5 × 5 grids until a result is found that matches the coordinates of the desired point within an acceptable tolerance. For the present work, the algorithm was required to determine the values of ξ and η that resulted in an error between the calculated coordinates and the actual coordinates of less than 1e-6. While there are more efficient numerical techniques for this part of the process, this approach converges reasonably quickly, usually within 6–8 iterations and is relatively inexpensive computationally. The approach can suffer some difficulty when the elements are significantly distorted. However, for the present work, the mesh was controlled sufficiently upfront and significant element distortions were avoided.

For a given point of interest in the grid that is being mapped to, once the parametric coordinates are known with an acceptable accuracy, any desirable field quantities can then be calculated. For this work, the field quantities of interest namely, displacements in the vertical and horizontal directions, the following relationships are used to compute those values:

$${(u}_{p};{v}_{p})=\sum\nolimits_{i=1}^{4}{N}_{i}\left(\xi ,\eta \right){(u}_{i};{v}_{i}).$$

(20)

This method is particularly advantageous because it avoids issues with averaging or smoothing around the crack tip or across the crack faces in the source data. This is because the target grid points are associated with elements from the output data. The nodal connectivity for the source elements is inherited from the source finite element model. Since the original mesh is created without elements spanning the crack tip or bridging across the crack faces, no averaging occurs due to target nodes on one side of the crack face being influenced by displacements of nodes on the opposing side of the crack. It should be noted that in the case of a set of points that are arranged in a rectangular fashion, this general method simplifies to bilinear interpolation.