- Email: [email protected]
A novel method for measuring discontinuous deformation in digital image correlation based on partition and dividing strategy
Accepted Manuscript A novel method for measuring discontinuous deformation in digital image correlation based on partition and dividing strategy Jing Han, Bing Pan PII: DOI: Reference:
S0013-7944(18)30467-3 https://doi.org/10.1016/j.engfracmech.2018.09.036 EFM 6171
To appear in:
Engineering Fracture Mechanics
Received Date: Revised Date: Accepted Date:
4 May 2018 28 September 2018 30 September 2018
Please cite this article as: Han, J., Pan, B., A novel method for measuring discontinuous deformation in digital image correlation based on partition and dividing strategy, Engineering Fracture Mechanics (2018), doi: https://doi.org/ 10.1016/j.engfracmech.2018.09.036
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A novel method for measuring discontinuous deformation in digital image correlation based on partition and dividing strategy Jing Han, Bing Pan* Institute of Solid Mechanics, Beihang University, Beijing 100191, China Email: [email protected] Abstract: Digital image correlation (DIC) is the most widely used photomechanics technique for surface deformation measurement. However, standard subset-based DIC restricts its application to continuous deformation, as the regularly-adopted shape functions fail to model discontinuous deformation with displacement jump when discontinuity presents. In this work, a simple but effective method based on the partition and division of interrogated subsets is presented to deal with discontinuity issues during DIC analysis. During the implementation of this method, correlation coefficient of each pixel within interrogated subsets is first calculated between the reference and deformed images to obtain a matrix of correlation coefficients. Then the correlation coefficient matrix is partitioned into two parts through binarization. Subsequently, a dividing line is fitted between the two sub-blocks, which can be considered as the initial discontinuity profile. After that, a differentiation strategy is adopted according to the relative sizes of the two sub-blocks, deformation parameters of which are retrieved using iterative Newton-Raphson method. The proposed method is verified using both numerically simulated images and real experimental images. Experimental results validate the efficacy and accuracy of the presented approach over existing methods in dealing with discontinuity issues. . Keywords: Digital image correlation, Discontinuous displacement fields, Fracture mechanics, Crack propagation
1. Introduction Digital image correlation (DIC) is a non-contact optical metrology for full-field shape and deformation measurement, which was originally proposed by Peters and Ranson in 1982 [1]. During the past three decades, DIC has undergone tremendous improvements and evolved to a versatile and powerful technique with prominent advantages of low-cost and simple setup, low requirements on experimental environments, and wide range of applicability [2-5]. In routine implementation of standard subset-based DIC, a region of interest (ROI) is first specified in the reference image, which is further divided into evenly spaced grid points. For each measurement point (i.e., grid point), a square subset of proper size centered at the considered point is chosen, and proper shape functions are defined to describe the underlying continuous deformation in deformed or target subsets. By optimizing the correlation function defined to quantify the similarity between the reference and the deformed subsets using a non-linear iterative optimization algorithm (e.g., Newton-Raphson algorithm), the displacement vector of each calculation point can be retrieved. Although the basic principle of standard subset-based DIC is straightforward and simple, close attention must be paid to the various details during its implementation. Presently, a zero-mean normalized sum of squared difference (ZNSSD), which is equivalent to the well-known zero-mean normalized cross-correlation (ZNCC) criterion, is recommended [6] to be used as the similarity measure in subset-based DIC due to its robustness to scale and offset changes of illumination as well as being easy for optimization. Specifically, a similarity measure can be established by combining ZNSSD criterion with proper shape functions as follows:
(1) where
is the desired deformation vector depicting the position and shape of the deformed targets to be sought. and denote gray scale levels in the reference image and deformed image, and , are the corresponding mean intensity values. For the regularly used first-order shape functions, which allow translation, rotation, shear, normal strains and their combinations of the subset, the deformation vector is written as , which can be determined by optimizing the non-linear correlation criterion defined in Eq.(1) using either classic Newton-Raphson algorithm [2] and the state-of-the-art inverse compositional Gauss-Newton (IC-GN) algorithm [7]. However, the abovementioned subset-based matching process faces challenges or even fails when dealing with specimens with discontinuity, such as cracks or shear bands [8-10]. This is because most of the existing subset-based DIC methods (even in commercial DIC software) inherently assumes that unknown displacement field is continuous and homogeneous in each subset, and the regularly-adopted first-order or second-order shape functions fail to model subsets that show strong discontinuities. As a result, the correlation quality dramatically deteriorates during image matching, leading to erroneous displacement and strain results. In earlier works to deal with discontinuities with DIC methods, displacements around discontinuities are either discarded due to low matching quality and inaccuracy or bypassed without even computing [11-14]. However, the displacement field around discontinuity is sometimes very valuable in the estimation of fracture mechanics parameters, for instance, it is essential to calculate stress intensity factor (SIF) or J integrals around the crack tip. Therefore, it is desirable to develop a simple yet effective DIC method to tackle discontinuity issues. In the literature, various improved DIC methods have been proposed to overcome this limitation [15-22]. For example, the subset-splitting method employed by Poissant et al [15] enables subsets to split into two sections along the crack. The two sections are then separately optimized to get in total twelve parameters, apart from two parameters modeling the crack as a line. Réthoré et al. introduced a finite element method into DIC (also known as extended-DIC), in which shape functions are enriched by Heaviside function to describe displacement jump [17, 18]. This method meshes the image with elements and computes the global correlation between reference and deformed images to get nodal displacements. Note that Heaviside function has also been introduced into local DIC shape functions to yield displacements around discontinuity. Recently, Hassan [21] proposed a discontinuous DIC (DDIC) method, which characterizes the discontinuity by direction of the tangent to the crack line and the corresponding Burgers vectors. In the subset-splitting method proposed by Poissant et al [15], the way to split a subset seems not stable and methods to obtain line parameters are not mentioned in detail. The accuracy of the splitting method adopted in this method relies heavily on initial guess. While the DDIC method takes a few times to iteratively optimize the angle of discontinuity line when there is discontinuity, which seems to be time-consuming and cumbersome. In short, the above-mentioned three methods adopt a similar strategy, though the global DIC differs fundamentally from local DIC, namely separate the subset (element) into two parts and model them with distinct parameters, which overcomes the disunity in deformation field caused by discontinuity. The challenging problem in above local discontinuous DIC methods involves the description and determination of discontinuity line. In addition to the above-mentioned three methods, there are also many alternative approaches adopting different strategies to tackle discontinuities. Nguyen et al [19] proposed a post-processing method, which recognizes the discontinuity point by a criterion based on the differences in displacements. Interpolation is then employed along the direction perpendicular to the discontinuity line. This method relies heavily on the accuracy of crack recognition. Jin et al [20] developed a pointwise DIC method, in which a genetic algorithm is employed to
cope with the large number of unknowns of in-plane displacements of all the pixels within a subset. However, all the methods mentioned above are either complicated, time-consuming or of low error tolerance. Here we propose a simple but effective method based on the partition and division of interrogated subsets to deal with discontinuity issues during DIC analysis. It combines the idea from the subset splitting method previously reported by Poissant et al [15] and logistic regression in machine learning [24]. As it takes full use of the surrounding information of each point in the computing subset, this method displays its superiority in recognizing the discontinuity line, which is more robust to noise than using matrix of square of displacement error in subset-splitting method [15]. Also, it enables line parameters to come out in a few iterations automatically without redundant computation. In addition, a differentiation strategy is adopted to reduce computation burden by optimizing parameters of only one sub-block when the two sub-blocks differ greatly in sizes. The proposed method is simple in principal, easy for implementation and can be considered as a supplement of the standard subset-based DIC. In the context below, we first elaborate the principles of proposed method, which is divided into three successive steps, namely subset partition and division as well as nonlinear optimization. To verify the effectiveness and accuracy of the presented approach over other methods in dealing with discontinuity issues, both numerically simulated images and real-world images were processed.
2. Principles Before introducing the proposed method, two assumptions are made first to simplify the discontinuous deformation problem: (1) discontinuity in the interrogated subset is recognized as a thin line, which is called discontinuity line in the following context; (2) the discontinuity line separates the interrogated subset into two sub-blocks that can be understood as different parts of the subset separated by discontinuity; and the corresponding two deformation fields are assumed to be different due to discontinuity. During the implementation of the proposed method, the interrogated subset is first divided along the discontinuity line into two sub-blocks, after which deformation vector is retrieved using iterative Newton-Raphson (NR) method. It should be noted that the proposed method is performed right after the implementation of standard subset-based DIC on those points, whose ZNCC correlation coefficients go below certain threshold (set to be0.9 in this work). In other words, recalculation is only required on the points with low correlation coefficients (as depicted in Fig. 1(a)). This section will introduce the method to partition and divide subsets where discontinuity crosses over, followed by the details of the refined optimization method. 2.1 Subset partition Partition is to identify crack locations in an interrogated subset and thus separates the subset into two parts (i.e., sub-blocks), while dividing is to figure out the parameters of crack line according to certain assumption. To partition an interrogated subset along discontinuity, a fact needs to be specified first that displacement jump due to discontinuity makes it impossible to transfer one accurate initial guess to both sub-blocks partitioned along discontinuity. In this case, only one sub-block can be well correlated where the point transferring initial guess locates, but the other sub-block could be poorly correlated. Based on this deduction, correlation coefficient of each pixel within the considered subset is calculated in the target subset using a small-sized subset (sub-window) centered at the interrogated pixel. To be more specific, let and denote gray scale value in the reference and deformed subsets respectively with , is initial guess obtained from neighboring point For any point in the considered subset, its counterpart in the deformed subset is . Thus correlation coefficient is calculated between and ( ) where sub-window size is pixels. As shown in Fig 1 (b), for an interrogated subset of pixels in recalculation area, correlation coefficients of a total of 961 pixels are calculated between and corresponding to obtain a
ZNCC matrix. For each pixel in the interrogated subset, a sub-window size of pixels is adopted for ZNCC calculation. Here the concept of sub-window size is used to distinguish it from subset size. Sub-window size is specially used for ZNCC matrix calculation, while subset size is used for correlation matching, though they are actually the same thing. As expected, correlation coefficients in the matrix will drop quickly where discontinuity presents, as shown in Fig 1(c). In order to extract this discontinuity line, the matrix of correlation coefficients is processed using binarization and divided into two sub-blocks along discontinuity. In this work, the two sub-blocks are called “zero-zone” and “one-zone”, respectively, where the binarization outputs are zeros and ones.
sub-window size
Calculation
subset size
Binarization
Fig.1 Illustration of the partition process: (a)
distribution within ROI after a pre-calculation using standard DIC; (b)
calculation for each point in the computing subset; (c)
distribution in the computing subset (d)
matrix
after binarization
2.2 Subset division Fig.1 (d) plots the ideal binary matrix, where the subset is completely partitioned into two sub-blocks. However, in most cases, “zero-zone” may include one values and “one-zone” may also include zero values. Even worse, the discontinuity line may not be very explicit due to the presence of noise. Based on these considerations, a
fitted dividing line is therefore more reasonable to cut the subset thoroughly. The dividing line can be assumed to be either straight or curved, depended on a trade-off between accuracy and efficiency. To be more explicit and without loss of generality, formulation of straight line and second-order curve are given as: (2) . (3) where are the parameters of the dividing line (three parameters are required to describe straight line whereas six parameters are needed for second-order curve); is the image coordinates in row and column direction, respectively. Although second-order curve fits the discontinuity line better, straight line assumption is preferred in this work for simplicity, as second-order approach seems to improve limited accuracy with remarkably increased computation time. To figure out the fitting parameters of the dividing line, logistic regression, which is a classic classification method in machine learning, is utilized [24]. As a regular linear model (Eq. (2)) outputs values greater than 1 or less than 0, running the linear function through a sigmoid activation function is a good choice to coerce reasonable answers. When a point is above the dividing line ( ), its real-valued input is mapped to 1 whereas when a point is below the dividing line ( ), its real-valued input is mapped to 0. A sigmoid function is given as: (4) Here
can be interpreted as the possibility when point
belongs to one-zone while
) denotes
the possibility when point belongs to zero-zone. is the true value of point in binary matrix of correlation coefficients. Then we can write the possibility of point in a compact form and get the objective function (loss function) . (5) (6) The loss function consists of two terms. When takes values of 0 or 1 for a given pixel, one of the two terms disappears. This loss function is commonly called log loss and is also commonly referred to as binary cross entropy. Various methods are available to maximize function to obtain parameters such as Newton-Raphson method or gradient ascent method. Here Newton-Raphson method is utilized to optimize objective function and iteration scheme is given in Eq. (8). =0 (7) (8) where denotes results obtained from a previous iteration, and initial guess set for Eq.(8) is given as [0, 0, 0 . First-order and second-order derivatives of the loss function L can be written in the following matrix forms. (9)
(10) 2.3 Displacement field determination around discontinuities After dividing the interrogated subset into two sub-blocks, deformation parameters at the subset center can be
determined using the classic iterative Newton-Raphson algorithm with incomplete subset size [23]. To this end, a differentiation strategy considering two cases is discussed here: Case #1: When one sub-block takes up little of the whole subset compared to the other, only the bigger sub-block is optimized. This is because pixel intensity values in the bigger sub-block suffices to yield accurate results, while the left smaller sub-block tends to reach an inferior solution for lack of information. ase #2: When there’s no significant difference in sizes, parameters of both sub-blocks are required to be optimized together, and deformation vector of the sub-block with higher correlation coefficient is taken as the final result. In this way, the effect of threshold taken in binarization procedure on dividing precision can be minimized. Also, the separate calculation of deformation vectors in sub-blocks avoids discontinuity problems arisen from cracks since deformation in each sub-block is continuous. However, to ensure accurate convergence of the Newton-Raphson algorithm for displacement estimation, proper initial guesses of the separated sub-blocks must be provided. Note that high-quality displacement field is available after standard DIC calculation except for area near discontinuity, thus initial guesses of the separated sub-blocks can be obtained from displacements at either side of the discontinuity. To initialize deformation vectors of both “one-zone” and “zero-zone” automatically and efficiently, the following strategy is adopted. Firstly correctly computed displacements of any neighboring point at any side of the discontinuity, whose ZNCC correlation coefficient is higher than 0.9, are first transferred to the interrogated point. Then this initial guess is assigned to the corresponding sub-block. Which sub-block a certain pixel belongs to is decided by Heaviside function with ， ， denoting line parameters obtained from logistic regression. Next, neighboring points of the interrogated point, whose displacements differ remarkably from the previous initial guess, are searched and transferred to the other sub-block. Searching range is limited here as failure points in standard DIC extends no more than four rows or columns in normal condition. (11) For easy explanation, the procedure for initial guess estimation is shown in Fig 2. Here is an interrogated point near crack. Displacements of point , whose ZNCC coefficient is larger than 0.9, is first transferred to “zero-zone” according to . Then by searching neighboring points of point
, point is found which satisfies and its displacement componentss are transferred to the other sub-block. Here the value of threshold is usually set as 1 pixel and can be adjusted according to actual displacement field. It should be noted that line parameters are solved before the implementation of Newton-Raphson method. However, if computed ZNCC correlation coefficients of the two sub-blocks still remain low (e.g., smaller than 0.85) after several iterations (e.g., 5 iterations), these line parameters will be updated by adjusting the threshold adopted in binarization of the correlation coefficient matrix in sub-block partition.
Crack
Fig 2. Process of estimating initial guess of two sub-blocks from displacement field obtained by standard DIC.
3. Verification by numerical simulation To examine the validity and performance of the proposed method, two sets of numerical simulation experiments were conducted. Computer generated speckle images with a resolution of 300×300 pixels are numerically deformed by two types of cracks including dislocation and semi-infinite type I crack both at different angles. 3.1 Simulation of dislocation Simulated speckle images are deformed by dislocation at an angle of and respectively, which is a constant discontinuity in displacement. As can be seen in Fig 3, horizontal dislocation originates from the center of the image and extends half of the image width until to the edge. Displacement field for horizontal dislocation is given as:
(12)
where is dislocation width, the Poisson’s ratio. Here Poisson’s ratio takes the value of 0.3 and which means displacement jump near discontinuity is 3 pixels and remains constant. (a)
Fig 3 Theoretical horizontal displacement field for dislocation at (a)
equals to 3,
(b)
and (b)
Standard DIC is performed first on the entire image and calculating points whose correlation coefficients go below 0.9 (referred to ) are considered near discontinuity as decorrelation appears in the region where displacement jump exists. With great possibility deformation vectors of these points are incorrect and need to be recalculated again using the proposed method. During the standard DIC implementation, optimization process stops once iteration times exceeds 20 or parameter increment . Displacements of 2704( ) points are calculated with a subset size of pixels and a 5-pixel grid step. After that, the proposed method is executed on the recalculation region. Generally, small subsets are preferred for discontinuous cases as deformation is usually complicated. Also, straight line assumption for dividing subsets into two sub-blocks requires smaller subset as discontinuity line is usually winding, in which case only when adopting small subsets curve can be reasonably regarded as straight. However, as interrogated subsets are partitioned along discontinuity, a subset size of pixels is adopted to guarantee sufficient intensity information. In addition, a sub-window size of pixels during the calculation of correlation coefficients matrix is used. Note that slightly lager or smaller sub-window sizes will not affect the performance greatly, but too large sub-window (larger than pixels) will probably blur discontinuity boundary while too small sub-window (smaller than pixels) lacks enough information to be robust to noise.
Figure 4 shows the theoretical and measured displacements of only ±1 pixel around dislocation. As shown in Fig 4 (c) and (f), correlation coefficients near discontinuity obtained by the standard DIC procedure drop dramatically, which is further reflected on the wrongly dispersed vertical and horizontal displacements (see blue up triangle) in Fig 4 (a) ,(b) and (d),(e). After a recalculation using the proposed method, the matching accuracy improves as subsets near dislocation are properly divided. As a result, displacement jump is well managed, as shown in Fig 4 (a) for horizontal dislocation and Fig 4 (d) for the inclined case. For the proposed method (see olive diamond), dislocation width roughly remains 3 pixels with 1.5-pixel displacement for the upper face and -1.5-pixel displacement for the lower face, which is in good agreement with theoretical results. (d)
1.5
standard DIC proposed method —— theoretical
1.0 0.5 0.0 -0.5 -1.0
1.5
U dispalcement (pixel)
U displacement (pixel)
(a)
1.0
standard DIC proposed method —— theoretical
0.5 0.0 -0.5 -1.0 -1.5
-1.5 0
50
100
150
200
250
300
0
50
100
X (pixel)
(e)
-0.4
standard DIC proposed method —— theoretical
-0.6 -0.8 -1.0 -1.2 -1.4 -1.6
V dispalcement (pixel)
V displacement (pixel)
(b)
150
200
250
300
X (pixel)
-0.4
standard DIC proposed method —— theoretical
-0.6 -0.8 -1.0 -1.2 -1.4 -1.6
-1.8
-1.8 0
50
100
150
200
250
300
0
50
100
X (pixel)
150
200
250
300
X (pixel)
(f) 1.0
(c) 1.0
ZNCC coefficient
ZNCC coefficient
0.9 0.9
0.8
0.8
0.7
standard DIC proposed method 0.7
standard DIC proposed method 0.6
50
100
150
200
250
300
0
50
X (pixel)
150
200
250
300
X (pixel)
Fig.4. Data collected along the dislocation line at ±1 pixels from dislocation at dislocation and (d), (e) and (f) are of
100
and
. Fig. (a), (b) and (c) are results of horizontal
dislocation, where (a) and (d) describe horizontal displacements, (b) and (e) vertical displacements, (c) and (f) ZNCC coefficients.
3.2 Simulation of mode I crack To further verify the effectiveness of the proposed method, mode I crack is then added to the reference image horizontally and at an angle of 45°respectively. Note that both the crack tips are located in the center of the image and crack extends from center to the edge of the image, whose span is half of the image width for horizontal crack and
fold of image width for 45°crack. The displacement fields are shown in Fig. 5 and theoretical equations of
horizontal mode I crack were given as:
(13)
where is the mode I stress intensity factor, is the Young’s modulus , is the Poisson’s ratio. ( ) is the polar coordinate originating from crack tip with denoting the distance from the crack tip, the angle relative to the horizontal axis. In order to clearly present the displacement fields, maximum displacement of roughly ，
pixels is guaranteed. And following values are taken:
，
.
(a)
(b)
Fig.5. Theoretical vertical displacement field for mode I crack at (a)
and (b)
.
By processing the simulated images using the proposed method, accuracy-enhanced displacement fields were retrieved. Fig 6 shows the horizontal and vertical displacements as well as the correlation coefficients of pixels near crack. Expected results are also obtained with the proposed method. As is seen in Fig 6, correlation coefficients (i.e., ) increase remarkably (basically larger than 0.9) compared to standard DIC, thus convincing displacements are retrieved with displacement jump in accordance with theoretical one. (d) 2.0 1.6
standard DIC proposed method —— theoretical
1.2
U displacement (pixel)
U displacement
(a)
0.8 0.4 0.0
proposed method standard DIC —— theoretical
1.5 1.0 0.5 0.0 -0.5
0
50
100
150
X (pixel)
200
250
300
0
50
100
150
X (pixel)
200
250
300
(e) 6 standard DIC proposed method —— theoretical
V displacement (pixel)
4 2 0 -2
standard DIC proposed method —— theoretical
4
V displacement (pixel)
(b)
2 0 -2 -4
-4 -6 0
50
100
150
200
250
300
0
50
100
X (pixel)
150
200
250
300
200
250
300
X (pixel)
(f) 1.0
(c) 1.0
0.8 0.7 0.6
standard DIC proposed method
ZNCC coefficient
ZNCC coefficient
0.9
0.9
0.8
0.7
proposed method standard DIC 0.6
0.5 0
50
100
150
200
250
0
300
Fig.6. Data collected along the crack line at ±1 pixels from crack at (e), (f) are of
50
100
150
X (pixlel)
X (pixel)
and
. Fig (a) (b) (c) are results of horizontal crack and (d),
crack where a, d describe horizontal displacements, (b), (e) vertical displacements and (c), (f) ZNCC coefficients.
3.3 Error analysis To quantitatively evaluate the calculated displacements, mean bias error and standard deviation error are adopted. Mean bias error of the measured displacement is defined as: (14) where , denotes the calculated and actual imposed sub-pixel displacement of the th point respectively. Standard deviation error of the measured displacement can be defined as: (15) Table 1 gives the mean bias error and standard deviation values around discontinuity for all the discontinuity types mentioned above, namely dislocation at and and semi-infinite mode I crack at and respectively. As shown in Table 1, after processing discontinuous area with the proposed method, both mean bias error and standard deviation drop greatly, which means results with better accuracy are obtained. In fact, on the whole, values of both mean bias and standard deviation will increase slightly around discontinuity compared to the entire image. This is because fewer intensity values of pixels are utilized during optimization, let alone partition bias when splitting subsets into two sub-blocks. Table 1 Mean bias error and standard deviation of displacement field for dislocation and mode I crack at areas around discontinuity
Standard DIC
Simulation result
The proposed method
dislocation at 0°
-0.3416 ± 1.2718
-0.0584 ± 0.0090
0.0875 ± 0.0503
0.0299 ± 0.0503
dislocation at 45º
0.0851 ± 0.4453
-0.0129 ± 0.0082
0.0192 ± 0.1133
0.0089 ± 0.0040
crack at 0°
-0.1067 ± 0.02737
0.2636 ± 1.9632
0.0467 ± 0.0432
-0.0391 ± 0.0766
crack at 45º
0.0380 ± 0.0358
0.3059 ± 2.2048
0.0214 ± 0.0263
-0.0294 ± 0.0417
4. Experimental verification To further demonstrate the validity of the proposed method, experiment was undertaken on fatigue cracks emanating from notches. Pure mode I loading condition was created ahead of a 23mm notch in a 5 mm thick CT specimen made of 7010 aluminum alloy. Specimen dimensions were in accordance to the guidelines found in the ASTM-E2472 and are illustrated in Fig 8(a) where unit is millimeter [30]. Loading frequency was at 15 Hz during the crack growth stages but the process was paused at intervals to allow for measurements to be taken at the various stages of crack growth. The recorded digital images are of 1024×1280 pixels in resolution and 256 gray levels in grayscale depth. During the calculation, a region of interest of 140×560 pixels is chosen and four consecutive images during loading are selected as the deformed images ,which are taken at 300min, 310min, 325 min, 350 min respectively at the load of 2100N. (a)
(b)
(c) ROI
(d)
Fig. 7 (a), (b), (c), (d) are images taken during the fracture experiment at 300 min, 310 min, 325 min and 350 min, respectively.
(a) (b)
Fig. 8 (a) The specimen configuration, (b) Force-time curve and red squares are when pictures are taken.
Totally, displacements of 2080(=20×104) points are calculated with a grid step of 5 pixels and a subset size of 21×21 pixels. Points whose correlation coefficients go below 0.95 (referred to ) are further processed with the proposed method with a 13×13 pixels sub-window size and a 31×31 pixels subset size. Data collected at ±2pixels along the crack line of the four deformed images is displayed in Fig 9. Through ZNCC coefficient results, we can see when crack opening displacement is relatively small, standard DIC can still achieve good results. With the growth of crack, standard DIC fails with ZNCC coefficients dropping remarkably. After processing with the proposed method, ZNCC coefficient improves greatly, which means displacement results near crack are reliable. As Newton-Raphson procedure probably don’t converge in standard DI , only vertical displacements obtained from the proposed method are shown in Fig 9. As we can see, displacement jump is successfully reconstructed by the proposed method, which becomes larger and larger with crack propagation. Also, crack tip position can be roughly estimated from the furcation of displacements near crack, which facilitate the accurate positioning of crack tip using the method proposed by Roux et all [29]. (b) 1.00
(a) 1.2 +2 pixel line +1 pixel line -1 pixel line -2 pixel line
0.8 0.6 0.4
0.95
ZNCC
V displacement (pixel)
1.0
x=254.56578
0.90
0.85
0.2 0.0
the proposed method standard DIC
0.80
-0.2 0
100
200
300
X (pixel)
400
500
0
100
200
300
X (pixel)
400
500
(d) 1.0 +2 pixel line +1 pixel line -1 pixel line -2 pixel line
1.5
0.9
1.0
ZNCC
V displacement (pixel)
(c) 2.0
0.5
0.8
x=311.40585
0.0
the proposed method standard DIC
0.7 -0.5 0
100
200
300
400
500
0
100
200
X (pixel)
+2 pixel line +1 pixel line -1 pixel line -2 pixel line
500
2.0 1.5 1.0
x=370.43207
0.9 0.8
ZNCC
2.5
0.7 0.6
the proposed method standard DIC
0.5
0.5
0.0 -0.5
0.4 0
100
200
300
400
0
500
100
200
-2 pixel line -1 pixel line +1 pixel line +2 pixel line
2
x=387.37747
1
0
0.8 0.7 0.6
standard DIC the proposed method
0.4 100
200
300
X (pixel)
500
0.9
0.5
0
400
(h) 1.0
ZNCC
(g) 4 3
300
X (pixel)
X (pixel)
V displacement (pixel)
400
(f) 1.0
(e) 3.0 V displacement (pixel)
300
X (pixel)
400
500
0
100
200
300
400
500
600
X (pixel)
Fig. 9 Data collected at ±2pixels along the crack line of images taken at 300min, 310min, 325 min, 350 min where (a), (b) are results of 300 min, (c), (d) are results of 310 min, (e), (f) are results of 325min and (g), (h) are results of 350min. Fig (a), (c), (e), (g) are vertical displacement results and (b), (d), (f), (h) are results of ZNCC coefficients.
Stress intensity factor (SIF) is then estimated to verify the correctness of the calculated displacement filed by the proposed method. William’s method is one of the most frequently used methods to estimate SIF from the displacement data of DIC [26-28].Through a least squares optimization, coefficients of Willian expansion terms can be determined, which have a direct relation with , and T stress. Besides, the crack tip position is determined by the ratio of two of the William coefficients mentioned in [29], which converges to an optimal solution. To validate the accuracy of the computing SIF value, theoretical solution is calculated in comparison with computing results. Theoretical solution is given in following equations where is crack length, is specimen thickness, P is loading force, is the distance from the hole center to specimen edge, which is 60 mm here.
(16) (17) (18) Results of SIF estimation and crack tip determination of the four deformed images are given in Table 2.As we can see, an error within the range of 5% is achieved through the displacement field obtained from the proposed method. Table 2 SIF estimation and crack tip determination of the four deformed images Image no.
Loading time /min
1 2 3 4
300 310 325 350
Initial crack tip location /pixel 232.70 311.41 370.43 387.37
Final crack tip location /pixel 246.7 311.1 384.7 413.4
Crack length /mm
Computing SIF /
11.18 13.98 17.17 18.42
9.78 11.25 14.28 15.90
m
Theoretical SIF /
m 10.16 11.13 14.80 15.69
Error /% 3.77 1.08 3.50 1.33
5. Conclusion A simple but effective method is proposed to deal with discontinuity problems in DIC. The proposed method can be divided into three successive steps, namely subset partition and division as well as nonlinear optimization. These three parts are first explained in detail. Simulation experiments were then performed to validate the performance of the proposed method by comparing the imposed displacement field with the computed one using images deformed by dislocation and semi-infinite mode I crack at the angle of 0º and 45º. Images from real experiments performed on an aluminum alloy are also analyzed. Results from both numerical simulation and real experiment show the excellent performance of the proposed method to reconstruct displacement field successfully around discontinuity. Based on the displacement field obtained near discontinuity, fracture mechanics parameters such as crack opening displacement, can be calculated directly. Furthermore, extracted from displacement field, other parameters such as stress intensity factor (SIF), crack tip position, crack length could be determined with higher accuracy, which are beneficial to quantitative crack growth monitoring.
Acknowledgements: This work is supported by the National Natural Science Foundation of China (Grant nos. 11872009 and 11632010), the National Key Research and Development Program of China (2018YFB0703500), the Aeronautical Science Foundation of China (2016ZD51034), and State Key Laboratory of Traction Power of Southwest Jiaotong University (Grant no.TPL1607).
Reference [1]. Peters WH, Ranson WF. Digital imaging techniques in experimental stress analysis. Opt. Eng.1981, 21:427-431. [2]. Bruck HA, McNeil SR, Sutton MA, Peters WH. Digital image correlation using Newton-Raphson method of partial differential correction. Experimental Mechanics.1989, 29(3): 261-267. [3]. Pan B, Qian KM, Xie HM, Asundi A. Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review. Measurement Science and Technology.2009, 20: 062001. [4]. Pan B. Digital image correlation for surface deformation measurements: historical developments, recent advances and future goals. Measurement Science and Technology.2018, 29(8): 082001. [5]. Pan B, Asundi A,Xie HM, Gao JX. Digital image correlation using iterative least squares and pointwise least squares for
displacement field and strain field measurements. Opt Lasers Eng.2008, 47(7-8):865-74 [6]. Pan B, Xie HM, Wang ZY. Equivalence of digital image correlation criteria. Applied Optics.2009, 49(28):5501-5509 [7]. Pan B, Li K, Tong W. Fast, robust and accurate digital image correlation calculation without redundant computations. Exp Mech. 2013, 53(7):1277-1289 [8]. McNeill SR, Peters WH, Sutton MA. Estimation of stress intensity factor by digital image correlation. Eng Fract Mech.1987, 281:101-112 [9]. Chasiotis I, Cho SW, Jonnalagadda K. Fracture toughness and subcritical crack growth in polycrystalline silicon. ASME.2006, 73:714–722 [10]. Rechenmacher AL. Grain-scale processes governing shear band initiation and evolution in sands. J Mech Phys Solids.2006,541:22–45 [11]. Sutton MA, Turner JL, Chao YJ, Bruck HA, Chae TL. Experimental investigations of three-dimensional effects near a crack tip using computer vision. Int J Fract. 1992,53(3):201–228 [12]. Gonzalez J, Knauss WG. Strain inhomogeneity and discontinuous crack growth in a particulate composite. J Mech Phys Solids. 1998,4610:1981–1995 [13]. Chao YJ, Luo PF, Kalthoff JF. An experimental study of the deformation fields around a propagating crack tip. Exp Mech. 1998,382:79–85 [14]. Barthelat F, Espinosa HD. An experimental investigation of deformation and fracture of nacre-mother of pearl. Exp Mech. 2007,473:311–324 [15]. Poissant J, Barthelat F. A novel subset splitting procedure for digital image correlation on discontinuous displacement fields. Exp Mech 2010, 50(3):353–64. [16]. V Valle, S Hedan, P. Cosenza, A.L. Fauchille, M. Berdjane. Digital Image Correlation Development for the Study of Materials Including Multiple Crossing Cracks. Exp Mech 2015, 55(2):379-391 [17]. Rthor J, Roux S, Hild F. From pictures to extended finite elements: extended digital image correlation (x-dic). CR Mcan 2007, 335(3):131-7. [18]. Rthor J, Tinnes J-P, Roux S, Buffire J-Y, Hild F. Extended three-dimensional digital image correlation (X3D-DIC). CR Mcan 2008, 336(8):643-9. [19]. Nguyen TL, Hall SA, Vacher P, Viggiani G. Fracture mechanisms in soft rock: identification and quantification of evolving displacement
discontinuities
by
extended
digital
image
correlation.
Tectonophysics
2011;503(12):117-28
[Thermo-Hydro-Chemo-Mechanical Couplings in Rock Physics and RockMechanic] [20]. Jin H, Bruck HA. Pointwise digital image correlation using genetic algorithms. Exp Tech. 2005,291:36-39 [21]. Ghulam Mubashar Hassan, Arcady V. Dyskin, Cara K. MacNish, Elena Pasternak, Igor Shufrin. Discontinuous Digital Image Correlation to reconstruct displacement and strain fields with discontinuities: Dislocation approach. Engineering Fracture Mechanics.2018,189:273-292 [22]. Helm J. Digital image correlation for specimens with multiple growing cracks. Exp Mech 2008, 48(6):753-62. [23]. Pan B, Wang ZY, Lu ZX. Genuine Full-field Deformation Measurement of An Object with Complex Shape using Reliability-guided Digital Image Correlation”. Optics Express. 2010, 18(2):1011-1023 [24]. https://en.wikipedia.org/wiki/Logistic_regression [25]. Mokhtarishirazabad M, Lopez-Crespo P, Moreno B, Lopez Moreno A, Zanganeh M. Evaluation of crack-tip fields from DIC data: A parametric study. Int J Fatigue. 2016, 89:11-19. [26]. Chiang FP, Asundi A. A white light speckle method applied to the determination of stress intensity factors and displacement field around a crack tip. Eng Fract Mech. 1981, 15(1-2):115‐ 121. [27]. Yoneyama S, Ogawa T, Kobayashi Y. Evaluating mixed-mode stress intensity factors from full-field displacement fields obtained by optical methods. Eng Fract Mech.2007, 74(9):1399‐ 1412. [28]. Yates JR, Zanganeh M, Tai YH. Quantifying crack tip displacement fields with DIC. Eng Fract Mech.2010,77(11): 2063-2076
[29]. Roux S, Réthoré J, Hild F. Digital image correlation and fracture: An advanced technique for estimating stress intensity factors of 2D and 3D cracks. J Phys D Appl Phys.2009, 42(21):214004. [30]. ASTM-E2472. Standard test method for determination of resistance to stable crack extension under low-constraint conditions; 2006.
Highlights
A simple but effective method based on the sub-block partition and area division strategy is presented to deal with discontinuity issues during DIC analysis.
The method provides a new method to partition the subsets near discontinuity into two parts, which is easy to implementation and robust to noise.
The method recognizes the discontinuity line using logistic regression and enables line parameters to come out in a few iterations automatically without redundant computation.
An improved shape function is adopted to describe displacement jump and strain variation in subsets near discontinuity to give a more precise description of discontinuous deformation.
S0013-7944(18)30467-3 https://doi.org/10.1016/j.engfracmech.2018.09.036 EFM 6171
To appear in:
Engineering Fracture Mechanics
Received Date: Revised Date: Accepted Date:
4 May 2018 28 September 2018 30 September 2018
Please cite this article as: Han, J., Pan, B., A novel method for measuring discontinuous deformation in digital image correlation based on partition and dividing strategy, Engineering Fracture Mechanics (2018), doi: https://doi.org/ 10.1016/j.engfracmech.2018.09.036
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A novel method for measuring discontinuous deformation in digital image correlation based on partition and dividing strategy Jing Han, Bing Pan* Institute of Solid Mechanics, Beihang University, Beijing 100191, China Email: [email protected] Abstract: Digital image correlation (DIC) is the most widely used photomechanics technique for surface deformation measurement. However, standard subset-based DIC restricts its application to continuous deformation, as the regularly-adopted shape functions fail to model discontinuous deformation with displacement jump when discontinuity presents. In this work, a simple but effective method based on the partition and division of interrogated subsets is presented to deal with discontinuity issues during DIC analysis. During the implementation of this method, correlation coefficient of each pixel within interrogated subsets is first calculated between the reference and deformed images to obtain a matrix of correlation coefficients. Then the correlation coefficient matrix is partitioned into two parts through binarization. Subsequently, a dividing line is fitted between the two sub-blocks, which can be considered as the initial discontinuity profile. After that, a differentiation strategy is adopted according to the relative sizes of the two sub-blocks, deformation parameters of which are retrieved using iterative Newton-Raphson method. The proposed method is verified using both numerically simulated images and real experimental images. Experimental results validate the efficacy and accuracy of the presented approach over existing methods in dealing with discontinuity issues. . Keywords: Digital image correlation, Discontinuous displacement fields, Fracture mechanics, Crack propagation
1. Introduction Digital image correlation (DIC) is a non-contact optical metrology for full-field shape and deformation measurement, which was originally proposed by Peters and Ranson in 1982 [1]. During the past three decades, DIC has undergone tremendous improvements and evolved to a versatile and powerful technique with prominent advantages of low-cost and simple setup, low requirements on experimental environments, and wide range of applicability [2-5]. In routine implementation of standard subset-based DIC, a region of interest (ROI) is first specified in the reference image, which is further divided into evenly spaced grid points. For each measurement point (i.e., grid point), a square subset of proper size centered at the considered point is chosen, and proper shape functions are defined to describe the underlying continuous deformation in deformed or target subsets. By optimizing the correlation function defined to quantify the similarity between the reference and the deformed subsets using a non-linear iterative optimization algorithm (e.g., Newton-Raphson algorithm), the displacement vector of each calculation point can be retrieved. Although the basic principle of standard subset-based DIC is straightforward and simple, close attention must be paid to the various details during its implementation. Presently, a zero-mean normalized sum of squared difference (ZNSSD), which is equivalent to the well-known zero-mean normalized cross-correlation (ZNCC) criterion, is recommended [6] to be used as the similarity measure in subset-based DIC due to its robustness to scale and offset changes of illumination as well as being easy for optimization. Specifically, a similarity measure can be established by combining ZNSSD criterion with proper shape functions as follows:
(1) where
is the desired deformation vector depicting the position and shape of the deformed targets to be sought. and denote gray scale levels in the reference image and deformed image, and , are the corresponding mean intensity values. For the regularly used first-order shape functions, which allow translation, rotation, shear, normal strains and their combinations of the subset, the deformation vector is written as , which can be determined by optimizing the non-linear correlation criterion defined in Eq.(1) using either classic Newton-Raphson algorithm [2] and the state-of-the-art inverse compositional Gauss-Newton (IC-GN) algorithm [7]. However, the abovementioned subset-based matching process faces challenges or even fails when dealing with specimens with discontinuity, such as cracks or shear bands [8-10]. This is because most of the existing subset-based DIC methods (even in commercial DIC software) inherently assumes that unknown displacement field is continuous and homogeneous in each subset, and the regularly-adopted first-order or second-order shape functions fail to model subsets that show strong discontinuities. As a result, the correlation quality dramatically deteriorates during image matching, leading to erroneous displacement and strain results. In earlier works to deal with discontinuities with DIC methods, displacements around discontinuities are either discarded due to low matching quality and inaccuracy or bypassed without even computing [11-14]. However, the displacement field around discontinuity is sometimes very valuable in the estimation of fracture mechanics parameters, for instance, it is essential to calculate stress intensity factor (SIF) or J integrals around the crack tip. Therefore, it is desirable to develop a simple yet effective DIC method to tackle discontinuity issues. In the literature, various improved DIC methods have been proposed to overcome this limitation [15-22]. For example, the subset-splitting method employed by Poissant et al [15] enables subsets to split into two sections along the crack. The two sections are then separately optimized to get in total twelve parameters, apart from two parameters modeling the crack as a line. Réthoré et al. introduced a finite element method into DIC (also known as extended-DIC), in which shape functions are enriched by Heaviside function to describe displacement jump [17, 18]. This method meshes the image with elements and computes the global correlation between reference and deformed images to get nodal displacements. Note that Heaviside function has also been introduced into local DIC shape functions to yield displacements around discontinuity. Recently, Hassan [21] proposed a discontinuous DIC (DDIC) method, which characterizes the discontinuity by direction of the tangent to the crack line and the corresponding Burgers vectors. In the subset-splitting method proposed by Poissant et al [15], the way to split a subset seems not stable and methods to obtain line parameters are not mentioned in detail. The accuracy of the splitting method adopted in this method relies heavily on initial guess. While the DDIC method takes a few times to iteratively optimize the angle of discontinuity line when there is discontinuity, which seems to be time-consuming and cumbersome. In short, the above-mentioned three methods adopt a similar strategy, though the global DIC differs fundamentally from local DIC, namely separate the subset (element) into two parts and model them with distinct parameters, which overcomes the disunity in deformation field caused by discontinuity. The challenging problem in above local discontinuous DIC methods involves the description and determination of discontinuity line. In addition to the above-mentioned three methods, there are also many alternative approaches adopting different strategies to tackle discontinuities. Nguyen et al [19] proposed a post-processing method, which recognizes the discontinuity point by a criterion based on the differences in displacements. Interpolation is then employed along the direction perpendicular to the discontinuity line. This method relies heavily on the accuracy of crack recognition. Jin et al [20] developed a pointwise DIC method, in which a genetic algorithm is employed to
cope with the large number of unknowns of in-plane displacements of all the pixels within a subset. However, all the methods mentioned above are either complicated, time-consuming or of low error tolerance. Here we propose a simple but effective method based on the partition and division of interrogated subsets to deal with discontinuity issues during DIC analysis. It combines the idea from the subset splitting method previously reported by Poissant et al [15] and logistic regression in machine learning [24]. As it takes full use of the surrounding information of each point in the computing subset, this method displays its superiority in recognizing the discontinuity line, which is more robust to noise than using matrix of square of displacement error in subset-splitting method [15]. Also, it enables line parameters to come out in a few iterations automatically without redundant computation. In addition, a differentiation strategy is adopted to reduce computation burden by optimizing parameters of only one sub-block when the two sub-blocks differ greatly in sizes. The proposed method is simple in principal, easy for implementation and can be considered as a supplement of the standard subset-based DIC. In the context below, we first elaborate the principles of proposed method, which is divided into three successive steps, namely subset partition and division as well as nonlinear optimization. To verify the effectiveness and accuracy of the presented approach over other methods in dealing with discontinuity issues, both numerically simulated images and real-world images were processed.
2. Principles Before introducing the proposed method, two assumptions are made first to simplify the discontinuous deformation problem: (1) discontinuity in the interrogated subset is recognized as a thin line, which is called discontinuity line in the following context; (2) the discontinuity line separates the interrogated subset into two sub-blocks that can be understood as different parts of the subset separated by discontinuity; and the corresponding two deformation fields are assumed to be different due to discontinuity. During the implementation of the proposed method, the interrogated subset is first divided along the discontinuity line into two sub-blocks, after which deformation vector is retrieved using iterative Newton-Raphson (NR) method. It should be noted that the proposed method is performed right after the implementation of standard subset-based DIC on those points, whose ZNCC correlation coefficients go below certain threshold (set to be0.9 in this work). In other words, recalculation is only required on the points with low correlation coefficients (as depicted in Fig. 1(a)). This section will introduce the method to partition and divide subsets where discontinuity crosses over, followed by the details of the refined optimization method. 2.1 Subset partition Partition is to identify crack locations in an interrogated subset and thus separates the subset into two parts (i.e., sub-blocks), while dividing is to figure out the parameters of crack line according to certain assumption. To partition an interrogated subset along discontinuity, a fact needs to be specified first that displacement jump due to discontinuity makes it impossible to transfer one accurate initial guess to both sub-blocks partitioned along discontinuity. In this case, only one sub-block can be well correlated where the point transferring initial guess locates, but the other sub-block could be poorly correlated. Based on this deduction, correlation coefficient of each pixel within the considered subset is calculated in the target subset using a small-sized subset (sub-window) centered at the interrogated pixel. To be more specific, let and denote gray scale value in the reference and deformed subsets respectively with , is initial guess obtained from neighboring point For any point in the considered subset, its counterpart in the deformed subset is . Thus correlation coefficient is calculated between and ( ) where sub-window size is pixels. As shown in Fig 1 (b), for an interrogated subset of pixels in recalculation area, correlation coefficients of a total of 961 pixels are calculated between and corresponding to obtain a
ZNCC matrix. For each pixel in the interrogated subset, a sub-window size of pixels is adopted for ZNCC calculation. Here the concept of sub-window size is used to distinguish it from subset size. Sub-window size is specially used for ZNCC matrix calculation, while subset size is used for correlation matching, though they are actually the same thing. As expected, correlation coefficients in the matrix will drop quickly where discontinuity presents, as shown in Fig 1(c). In order to extract this discontinuity line, the matrix of correlation coefficients is processed using binarization and divided into two sub-blocks along discontinuity. In this work, the two sub-blocks are called “zero-zone” and “one-zone”, respectively, where the binarization outputs are zeros and ones.
sub-window size
Calculation
subset size
Binarization
Fig.1 Illustration of the partition process: (a)
distribution within ROI after a pre-calculation using standard DIC; (b)
calculation for each point in the computing subset; (c)
distribution in the computing subset (d)
matrix
after binarization
2.2 Subset division Fig.1 (d) plots the ideal binary matrix, where the subset is completely partitioned into two sub-blocks. However, in most cases, “zero-zone” may include one values and “one-zone” may also include zero values. Even worse, the discontinuity line may not be very explicit due to the presence of noise. Based on these considerations, a
fitted dividing line is therefore more reasonable to cut the subset thoroughly. The dividing line can be assumed to be either straight or curved, depended on a trade-off between accuracy and efficiency. To be more explicit and without loss of generality, formulation of straight line and second-order curve are given as: (2) . (3) where are the parameters of the dividing line (three parameters are required to describe straight line whereas six parameters are needed for second-order curve); is the image coordinates in row and column direction, respectively. Although second-order curve fits the discontinuity line better, straight line assumption is preferred in this work for simplicity, as second-order approach seems to improve limited accuracy with remarkably increased computation time. To figure out the fitting parameters of the dividing line, logistic regression, which is a classic classification method in machine learning, is utilized [24]. As a regular linear model (Eq. (2)) outputs values greater than 1 or less than 0, running the linear function through a sigmoid activation function is a good choice to coerce reasonable answers. When a point is above the dividing line ( ), its real-valued input is mapped to 1 whereas when a point is below the dividing line ( ), its real-valued input is mapped to 0. A sigmoid function is given as: (4) Here
can be interpreted as the possibility when point
belongs to one-zone while
) denotes
the possibility when point belongs to zero-zone. is the true value of point in binary matrix of correlation coefficients. Then we can write the possibility of point in a compact form and get the objective function (loss function) . (5) (6) The loss function consists of two terms. When takes values of 0 or 1 for a given pixel, one of the two terms disappears. This loss function is commonly called log loss and is also commonly referred to as binary cross entropy. Various methods are available to maximize function to obtain parameters such as Newton-Raphson method or gradient ascent method. Here Newton-Raphson method is utilized to optimize objective function and iteration scheme is given in Eq. (8). =0 (7) (8) where denotes results obtained from a previous iteration, and initial guess set for Eq.(8) is given as [0, 0, 0 . First-order and second-order derivatives of the loss function L can be written in the following matrix forms. (9)
(10) 2.3 Displacement field determination around discontinuities After dividing the interrogated subset into two sub-blocks, deformation parameters at the subset center can be
determined using the classic iterative Newton-Raphson algorithm with incomplete subset size [23]. To this end, a differentiation strategy considering two cases is discussed here: Case #1: When one sub-block takes up little of the whole subset compared to the other, only the bigger sub-block is optimized. This is because pixel intensity values in the bigger sub-block suffices to yield accurate results, while the left smaller sub-block tends to reach an inferior solution for lack of information. ase #2: When there’s no significant difference in sizes, parameters of both sub-blocks are required to be optimized together, and deformation vector of the sub-block with higher correlation coefficient is taken as the final result. In this way, the effect of threshold taken in binarization procedure on dividing precision can be minimized. Also, the separate calculation of deformation vectors in sub-blocks avoids discontinuity problems arisen from cracks since deformation in each sub-block is continuous. However, to ensure accurate convergence of the Newton-Raphson algorithm for displacement estimation, proper initial guesses of the separated sub-blocks must be provided. Note that high-quality displacement field is available after standard DIC calculation except for area near discontinuity, thus initial guesses of the separated sub-blocks can be obtained from displacements at either side of the discontinuity. To initialize deformation vectors of both “one-zone” and “zero-zone” automatically and efficiently, the following strategy is adopted. Firstly correctly computed displacements of any neighboring point at any side of the discontinuity, whose ZNCC correlation coefficient is higher than 0.9, are first transferred to the interrogated point. Then this initial guess is assigned to the corresponding sub-block. Which sub-block a certain pixel belongs to is decided by Heaviside function with ， ， denoting line parameters obtained from logistic regression. Next, neighboring points of the interrogated point, whose displacements differ remarkably from the previous initial guess, are searched and transferred to the other sub-block. Searching range is limited here as failure points in standard DIC extends no more than four rows or columns in normal condition. (11) For easy explanation, the procedure for initial guess estimation is shown in Fig 2. Here is an interrogated point near crack. Displacements of point , whose ZNCC coefficient is larger than 0.9, is first transferred to “zero-zone” according to . Then by searching neighboring points of point
, point is found which satisfies and its displacement componentss are transferred to the other sub-block. Here the value of threshold is usually set as 1 pixel and can be adjusted according to actual displacement field. It should be noted that line parameters are solved before the implementation of Newton-Raphson method. However, if computed ZNCC correlation coefficients of the two sub-blocks still remain low (e.g., smaller than 0.85) after several iterations (e.g., 5 iterations), these line parameters will be updated by adjusting the threshold adopted in binarization of the correlation coefficient matrix in sub-block partition.
Crack
Fig 2. Process of estimating initial guess of two sub-blocks from displacement field obtained by standard DIC.
3. Verification by numerical simulation To examine the validity and performance of the proposed method, two sets of numerical simulation experiments were conducted. Computer generated speckle images with a resolution of 300×300 pixels are numerically deformed by two types of cracks including dislocation and semi-infinite type I crack both at different angles. 3.1 Simulation of dislocation Simulated speckle images are deformed by dislocation at an angle of and respectively, which is a constant discontinuity in displacement. As can be seen in Fig 3, horizontal dislocation originates from the center of the image and extends half of the image width until to the edge. Displacement field for horizontal dislocation is given as:
(12)
where is dislocation width, the Poisson’s ratio. Here Poisson’s ratio takes the value of 0.3 and which means displacement jump near discontinuity is 3 pixels and remains constant. (a)
Fig 3 Theoretical horizontal displacement field for dislocation at (a)
equals to 3,
(b)
and (b)
Standard DIC is performed first on the entire image and calculating points whose correlation coefficients go below 0.9 (referred to ) are considered near discontinuity as decorrelation appears in the region where displacement jump exists. With great possibility deformation vectors of these points are incorrect and need to be recalculated again using the proposed method. During the standard DIC implementation, optimization process stops once iteration times exceeds 20 or parameter increment . Displacements of 2704( ) points are calculated with a subset size of pixels and a 5-pixel grid step. After that, the proposed method is executed on the recalculation region. Generally, small subsets are preferred for discontinuous cases as deformation is usually complicated. Also, straight line assumption for dividing subsets into two sub-blocks requires smaller subset as discontinuity line is usually winding, in which case only when adopting small subsets curve can be reasonably regarded as straight. However, as interrogated subsets are partitioned along discontinuity, a subset size of pixels is adopted to guarantee sufficient intensity information. In addition, a sub-window size of pixels during the calculation of correlation coefficients matrix is used. Note that slightly lager or smaller sub-window sizes will not affect the performance greatly, but too large sub-window (larger than pixels) will probably blur discontinuity boundary while too small sub-window (smaller than pixels) lacks enough information to be robust to noise.
Figure 4 shows the theoretical and measured displacements of only ±1 pixel around dislocation. As shown in Fig 4 (c) and (f), correlation coefficients near discontinuity obtained by the standard DIC procedure drop dramatically, which is further reflected on the wrongly dispersed vertical and horizontal displacements (see blue up triangle) in Fig 4 (a) ,(b) and (d),(e). After a recalculation using the proposed method, the matching accuracy improves as subsets near dislocation are properly divided. As a result, displacement jump is well managed, as shown in Fig 4 (a) for horizontal dislocation and Fig 4 (d) for the inclined case. For the proposed method (see olive diamond), dislocation width roughly remains 3 pixels with 1.5-pixel displacement for the upper face and -1.5-pixel displacement for the lower face, which is in good agreement with theoretical results. (d)
1.5
standard DIC proposed method —— theoretical
1.0 0.5 0.0 -0.5 -1.0
1.5
U dispalcement (pixel)
U displacement (pixel)
(a)
1.0
standard DIC proposed method —— theoretical
0.5 0.0 -0.5 -1.0 -1.5
-1.5 0
50
100
150
200
250
300
0
50
100
X (pixel)
(e)
-0.4
standard DIC proposed method —— theoretical
-0.6 -0.8 -1.0 -1.2 -1.4 -1.6
V dispalcement (pixel)
V displacement (pixel)
(b)
150
200
250
300
X (pixel)
-0.4
standard DIC proposed method —— theoretical
-0.6 -0.8 -1.0 -1.2 -1.4 -1.6
-1.8
-1.8 0
50
100
150
200
250
300
0
50
100
X (pixel)
150
200
250
300
X (pixel)
(f) 1.0
(c) 1.0
ZNCC coefficient
ZNCC coefficient
0.9 0.9
0.8
0.8
0.7
standard DIC proposed method 0.7
standard DIC proposed method 0.6
50
100
150
200
250
300
0
50
X (pixel)
150
200
250
300
X (pixel)
Fig.4. Data collected along the dislocation line at ±1 pixels from dislocation at dislocation and (d), (e) and (f) are of
100
and
. Fig. (a), (b) and (c) are results of horizontal
dislocation, where (a) and (d) describe horizontal displacements, (b) and (e) vertical displacements, (c) and (f) ZNCC coefficients.
3.2 Simulation of mode I crack To further verify the effectiveness of the proposed method, mode I crack is then added to the reference image horizontally and at an angle of 45°respectively. Note that both the crack tips are located in the center of the image and crack extends from center to the edge of the image, whose span is half of the image width for horizontal crack and
fold of image width for 45°crack. The displacement fields are shown in Fig. 5 and theoretical equations of
horizontal mode I crack were given as:
(13)
where is the mode I stress intensity factor, is the Young’s modulus , is the Poisson’s ratio. ( ) is the polar coordinate originating from crack tip with denoting the distance from the crack tip, the angle relative to the horizontal axis. In order to clearly present the displacement fields, maximum displacement of roughly ，
pixels is guaranteed. And following values are taken:
，
.
(a)
(b)
Fig.5. Theoretical vertical displacement field for mode I crack at (a)
and (b)
.
By processing the simulated images using the proposed method, accuracy-enhanced displacement fields were retrieved. Fig 6 shows the horizontal and vertical displacements as well as the correlation coefficients of pixels near crack. Expected results are also obtained with the proposed method. As is seen in Fig 6, correlation coefficients (i.e., ) increase remarkably (basically larger than 0.9) compared to standard DIC, thus convincing displacements are retrieved with displacement jump in accordance with theoretical one. (d) 2.0 1.6
standard DIC proposed method —— theoretical
1.2
U displacement (pixel)
U displacement
(a)
0.8 0.4 0.0
proposed method standard DIC —— theoretical
1.5 1.0 0.5 0.0 -0.5
0
50
100
150
X (pixel)
200
250
300
0
50
100
150
X (pixel)
200
250
300
(e) 6 standard DIC proposed method —— theoretical
V displacement (pixel)
4 2 0 -2
standard DIC proposed method —— theoretical
4
V displacement (pixel)
(b)
2 0 -2 -4
-4 -6 0
50
100
150
200
250
300
0
50
100
X (pixel)
150
200
250
300
200
250
300
X (pixel)
(f) 1.0
(c) 1.0
0.8 0.7 0.6
standard DIC proposed method
ZNCC coefficient
ZNCC coefficient
0.9
0.9
0.8
0.7
proposed method standard DIC 0.6
0.5 0
50
100
150
200
250
0
300
Fig.6. Data collected along the crack line at ±1 pixels from crack at (e), (f) are of
50
100
150
X (pixlel)
X (pixel)
and
. Fig (a) (b) (c) are results of horizontal crack and (d),
crack where a, d describe horizontal displacements, (b), (e) vertical displacements and (c), (f) ZNCC coefficients.
3.3 Error analysis To quantitatively evaluate the calculated displacements, mean bias error and standard deviation error are adopted. Mean bias error of the measured displacement is defined as: (14) where , denotes the calculated and actual imposed sub-pixel displacement of the th point respectively. Standard deviation error of the measured displacement can be defined as: (15) Table 1 gives the mean bias error and standard deviation values around discontinuity for all the discontinuity types mentioned above, namely dislocation at and and semi-infinite mode I crack at and respectively. As shown in Table 1, after processing discontinuous area with the proposed method, both mean bias error and standard deviation drop greatly, which means results with better accuracy are obtained. In fact, on the whole, values of both mean bias and standard deviation will increase slightly around discontinuity compared to the entire image. This is because fewer intensity values of pixels are utilized during optimization, let alone partition bias when splitting subsets into two sub-blocks. Table 1 Mean bias error and standard deviation of displacement field for dislocation and mode I crack at areas around discontinuity
Standard DIC
Simulation result
The proposed method
dislocation at 0°
-0.3416 ± 1.2718
-0.0584 ± 0.0090
0.0875 ± 0.0503
0.0299 ± 0.0503
dislocation at 45º
0.0851 ± 0.4453
-0.0129 ± 0.0082
0.0192 ± 0.1133
0.0089 ± 0.0040
crack at 0°
-0.1067 ± 0.02737
0.2636 ± 1.9632
0.0467 ± 0.0432
-0.0391 ± 0.0766
crack at 45º
0.0380 ± 0.0358
0.3059 ± 2.2048
0.0214 ± 0.0263
-0.0294 ± 0.0417
4. Experimental verification To further demonstrate the validity of the proposed method, experiment was undertaken on fatigue cracks emanating from notches. Pure mode I loading condition was created ahead of a 23mm notch in a 5 mm thick CT specimen made of 7010 aluminum alloy. Specimen dimensions were in accordance to the guidelines found in the ASTM-E2472 and are illustrated in Fig 8(a) where unit is millimeter [30]. Loading frequency was at 15 Hz during the crack growth stages but the process was paused at intervals to allow for measurements to be taken at the various stages of crack growth. The recorded digital images are of 1024×1280 pixels in resolution and 256 gray levels in grayscale depth. During the calculation, a region of interest of 140×560 pixels is chosen and four consecutive images during loading are selected as the deformed images ,which are taken at 300min, 310min, 325 min, 350 min respectively at the load of 2100N. (a)
(b)
(c) ROI
(d)
Fig. 7 (a), (b), (c), (d) are images taken during the fracture experiment at 300 min, 310 min, 325 min and 350 min, respectively.
(a) (b)
Fig. 8 (a) The specimen configuration, (b) Force-time curve and red squares are when pictures are taken.
Totally, displacements of 2080(=20×104) points are calculated with a grid step of 5 pixels and a subset size of 21×21 pixels. Points whose correlation coefficients go below 0.95 (referred to ) are further processed with the proposed method with a 13×13 pixels sub-window size and a 31×31 pixels subset size. Data collected at ±2pixels along the crack line of the four deformed images is displayed in Fig 9. Through ZNCC coefficient results, we can see when crack opening displacement is relatively small, standard DIC can still achieve good results. With the growth of crack, standard DIC fails with ZNCC coefficients dropping remarkably. After processing with the proposed method, ZNCC coefficient improves greatly, which means displacement results near crack are reliable. As Newton-Raphson procedure probably don’t converge in standard DI , only vertical displacements obtained from the proposed method are shown in Fig 9. As we can see, displacement jump is successfully reconstructed by the proposed method, which becomes larger and larger with crack propagation. Also, crack tip position can be roughly estimated from the furcation of displacements near crack, which facilitate the accurate positioning of crack tip using the method proposed by Roux et all [29]. (b) 1.00
(a) 1.2 +2 pixel line +1 pixel line -1 pixel line -2 pixel line
0.8 0.6 0.4
0.95
ZNCC
V displacement (pixel)
1.0
x=254.56578
0.90
0.85
0.2 0.0
the proposed method standard DIC
0.80
-0.2 0
100
200
300
X (pixel)
400
500
0
100
200
300
X (pixel)
400
500
(d) 1.0 +2 pixel line +1 pixel line -1 pixel line -2 pixel line
1.5
0.9
1.0
ZNCC
V displacement (pixel)
(c) 2.0
0.5
0.8
x=311.40585
0.0
the proposed method standard DIC
0.7 -0.5 0
100
200
300
400
500
0
100
200
X (pixel)
+2 pixel line +1 pixel line -1 pixel line -2 pixel line
500
2.0 1.5 1.0
x=370.43207
0.9 0.8
ZNCC
2.5
0.7 0.6
the proposed method standard DIC
0.5
0.5
0.0 -0.5
0.4 0
100
200
300
400
0
500
100
200
-2 pixel line -1 pixel line +1 pixel line +2 pixel line
2
x=387.37747
1
0
0.8 0.7 0.6
standard DIC the proposed method
0.4 100
200
300
X (pixel)
500
0.9
0.5
0
400
(h) 1.0
ZNCC
(g) 4 3
300
X (pixel)
X (pixel)
V displacement (pixel)
400
(f) 1.0
(e) 3.0 V displacement (pixel)
300
X (pixel)
400
500
0
100
200
300
400
500
600
X (pixel)
Fig. 9 Data collected at ±2pixels along the crack line of images taken at 300min, 310min, 325 min, 350 min where (a), (b) are results of 300 min, (c), (d) are results of 310 min, (e), (f) are results of 325min and (g), (h) are results of 350min. Fig (a), (c), (e), (g) are vertical displacement results and (b), (d), (f), (h) are results of ZNCC coefficients.
Stress intensity factor (SIF) is then estimated to verify the correctness of the calculated displacement filed by the proposed method. William’s method is one of the most frequently used methods to estimate SIF from the displacement data of DIC [26-28].Through a least squares optimization, coefficients of Willian expansion terms can be determined, which have a direct relation with , and T stress. Besides, the crack tip position is determined by the ratio of two of the William coefficients mentioned in [29], which converges to an optimal solution. To validate the accuracy of the computing SIF value, theoretical solution is calculated in comparison with computing results. Theoretical solution is given in following equations where is crack length, is specimen thickness, P is loading force, is the distance from the hole center to specimen edge, which is 60 mm here.
(16) (17) (18) Results of SIF estimation and crack tip determination of the four deformed images are given in Table 2.As we can see, an error within the range of 5% is achieved through the displacement field obtained from the proposed method. Table 2 SIF estimation and crack tip determination of the four deformed images Image no.
Loading time /min
1 2 3 4
300 310 325 350
Initial crack tip location /pixel 232.70 311.41 370.43 387.37
Final crack tip location /pixel 246.7 311.1 384.7 413.4
Crack length /mm
Computing SIF /
11.18 13.98 17.17 18.42
9.78 11.25 14.28 15.90
m
Theoretical SIF /
m 10.16 11.13 14.80 15.69
Error /% 3.77 1.08 3.50 1.33
5. Conclusion A simple but effective method is proposed to deal with discontinuity problems in DIC. The proposed method can be divided into three successive steps, namely subset partition and division as well as nonlinear optimization. These three parts are first explained in detail. Simulation experiments were then performed to validate the performance of the proposed method by comparing the imposed displacement field with the computed one using images deformed by dislocation and semi-infinite mode I crack at the angle of 0º and 45º. Images from real experiments performed on an aluminum alloy are also analyzed. Results from both numerical simulation and real experiment show the excellent performance of the proposed method to reconstruct displacement field successfully around discontinuity. Based on the displacement field obtained near discontinuity, fracture mechanics parameters such as crack opening displacement, can be calculated directly. Furthermore, extracted from displacement field, other parameters such as stress intensity factor (SIF), crack tip position, crack length could be determined with higher accuracy, which are beneficial to quantitative crack growth monitoring.
Acknowledgements: This work is supported by the National Natural Science Foundation of China (Grant nos. 11872009 and 11632010), the National Key Research and Development Program of China (2018YFB0703500), the Aeronautical Science Foundation of China (2016ZD51034), and State Key Laboratory of Traction Power of Southwest Jiaotong University (Grant no.TPL1607).
Reference [1]. Peters WH, Ranson WF. Digital imaging techniques in experimental stress analysis. Opt. Eng.1981, 21:427-431. [2]. Bruck HA, McNeil SR, Sutton MA, Peters WH. Digital image correlation using Newton-Raphson method of partial differential correction. Experimental Mechanics.1989, 29(3): 261-267. [3]. Pan B, Qian KM, Xie HM, Asundi A. Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review. Measurement Science and Technology.2009, 20: 062001. [4]. Pan B. Digital image correlation for surface deformation measurements: historical developments, recent advances and future goals. Measurement Science and Technology.2018, 29(8): 082001. [5]. Pan B, Asundi A,Xie HM, Gao JX. Digital image correlation using iterative least squares and pointwise least squares for
displacement field and strain field measurements. Opt Lasers Eng.2008, 47(7-8):865-74 [6]. Pan B, Xie HM, Wang ZY. Equivalence of digital image correlation criteria. Applied Optics.2009, 49(28):5501-5509 [7]. Pan B, Li K, Tong W. Fast, robust and accurate digital image correlation calculation without redundant computations. Exp Mech. 2013, 53(7):1277-1289 [8]. McNeill SR, Peters WH, Sutton MA. Estimation of stress intensity factor by digital image correlation. Eng Fract Mech.1987, 281:101-112 [9]. Chasiotis I, Cho SW, Jonnalagadda K. Fracture toughness and subcritical crack growth in polycrystalline silicon. ASME.2006, 73:714–722 [10]. Rechenmacher AL. Grain-scale processes governing shear band initiation and evolution in sands. J Mech Phys Solids.2006,541:22–45 [11]. Sutton MA, Turner JL, Chao YJ, Bruck HA, Chae TL. Experimental investigations of three-dimensional effects near a crack tip using computer vision. Int J Fract. 1992,53(3):201–228 [12]. Gonzalez J, Knauss WG. Strain inhomogeneity and discontinuous crack growth in a particulate composite. J Mech Phys Solids. 1998,4610:1981–1995 [13]. Chao YJ, Luo PF, Kalthoff JF. An experimental study of the deformation fields around a propagating crack tip. Exp Mech. 1998,382:79–85 [14]. Barthelat F, Espinosa HD. An experimental investigation of deformation and fracture of nacre-mother of pearl. Exp Mech. 2007,473:311–324 [15]. Poissant J, Barthelat F. A novel subset splitting procedure for digital image correlation on discontinuous displacement fields. Exp Mech 2010, 50(3):353–64. [16]. V Valle, S Hedan, P. Cosenza, A.L. Fauchille, M. Berdjane. Digital Image Correlation Development for the Study of Materials Including Multiple Crossing Cracks. Exp Mech 2015, 55(2):379-391 [17]. Rthor J, Roux S, Hild F. From pictures to extended finite elements: extended digital image correlation (x-dic). CR Mcan 2007, 335(3):131-7. [18]. Rthor J, Tinnes J-P, Roux S, Buffire J-Y, Hild F. Extended three-dimensional digital image correlation (X3D-DIC). CR Mcan 2008, 336(8):643-9. [19]. Nguyen TL, Hall SA, Vacher P, Viggiani G. Fracture mechanisms in soft rock: identification and quantification of evolving displacement
discontinuities
by
extended
digital
image
correlation.
Tectonophysics
2011;503(12):117-28
[Thermo-Hydro-Chemo-Mechanical Couplings in Rock Physics and RockMechanic] [20]. Jin H, Bruck HA. Pointwise digital image correlation using genetic algorithms. Exp Tech. 2005,291:36-39 [21]. Ghulam Mubashar Hassan, Arcady V. Dyskin, Cara K. MacNish, Elena Pasternak, Igor Shufrin. Discontinuous Digital Image Correlation to reconstruct displacement and strain fields with discontinuities: Dislocation approach. Engineering Fracture Mechanics.2018,189:273-292 [22]. Helm J. Digital image correlation for specimens with multiple growing cracks. Exp Mech 2008, 48(6):753-62. [23]. Pan B, Wang ZY, Lu ZX. Genuine Full-field Deformation Measurement of An Object with Complex Shape using Reliability-guided Digital Image Correlation”. Optics Express. 2010, 18(2):1011-1023 [24]. https://en.wikipedia.org/wiki/Logistic_regression [25]. Mokhtarishirazabad M, Lopez-Crespo P, Moreno B, Lopez Moreno A, Zanganeh M. Evaluation of crack-tip fields from DIC data: A parametric study. Int J Fatigue. 2016, 89:11-19. [26]. Chiang FP, Asundi A. A white light speckle method applied to the determination of stress intensity factors and displacement field around a crack tip. Eng Fract Mech. 1981, 15(1-2):115‐ 121. [27]. Yoneyama S, Ogawa T, Kobayashi Y. Evaluating mixed-mode stress intensity factors from full-field displacement fields obtained by optical methods. Eng Fract Mech.2007, 74(9):1399‐ 1412. [28]. Yates JR, Zanganeh M, Tai YH. Quantifying crack tip displacement fields with DIC. Eng Fract Mech.2010,77(11): 2063-2076
[29]. Roux S, Réthoré J, Hild F. Digital image correlation and fracture: An advanced technique for estimating stress intensity factors of 2D and 3D cracks. J Phys D Appl Phys.2009, 42(21):214004. [30]. ASTM-E2472. Standard test method for determination of resistance to stable crack extension under low-constraint conditions; 2006.
Highlights
A simple but effective method based on the sub-block partition and area division strategy is presented to deal with discontinuity issues during DIC analysis.
The method provides a new method to partition the subsets near discontinuity into two parts, which is easy to implementation and robust to noise.
The method recognizes the discontinuity line using logistic regression and enables line parameters to come out in a few iterations automatically without redundant computation.
An improved shape function is adopted to describe displacement jump and strain variation in subsets near discontinuity to give a more precise description of discontinuous deformation.