- Email: [email protected]
A weak-form method for strain field calculation with digital image correlation
Optics and Lasers in Engineering 121 (2019) 495–505
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
A weak-form method for strain field calculation with digital image correlation Leiting Dong a,∗, Qiuhong Yang b, Junbo Wang a a b
School of Aeronautic Science and Engineering, Beihang University, Beijing 100191, China COMAC Shanghai Aircraft Design and Research Institute, Shanghai 201210, China
a r t i c l e
i n f o
Keywords: Digital image correlation Strain field measurement Denoising Weak-form Test functions Trial functions
a b s t r a c t In digital image correlation (DIC), iterative spatial-domain cross-correlation algorithms have been routinely used to extract displacement fields from the recorded images. Then, strain fields are computed from the noisy displacement fields by a proper numerical different approach, like pointwise least squares fitting, which is tricky in choosing optimal calculation parameters. This work proposes an alternative weak-form framework to calculate strain fields from noisy displacement measurements by DIC. While direct numerical differentiation may amplify the noise in the displacement fields, the proposed method transfers differentiation into integration by employing the idea of the weak-form, which is the foundation of modern computational mechanics. By selecting cosine series as trial/test functions and employing their orthogonality, a semi-analytical formula for calculating the strain field is given. Moreover, the Filon’s method is used in this study to improve the accuracy of numerical integration as sinusoidal terms are involved. The effectiveness of the proposed method is verified by numerical examples. It is found that the proposed method is comparable to the famous PLS algorithm when it is for simple strain fields, while it is better in accuracy and its insensitivity to noise when large-gradient or highly-oscillating strain fields are to be measured.
1. Introduction Surface deformation measurement of materials and structures is an important task of experimental mechanics. Various optical methods, such as holography interferometry, electronic speckle pattern interferometry (ESPI) [1], digital image correlation (DIC) [2,3], have been proposed and advocated for non-contact and full-field deformation measurement. Among these methods, DIC techniques proposed by Peters et al. [4,5] have been extensively investigated and used due to its simplicity, robustness and wide range of applicability [2]. The basic principle of DIC method is simple and straightforward, namely matching the subsets of the surface in the reference image and the deformed image to calculate the displacement fields and strain fields. To this purpose, certain correlation criteria [2], e.g., the zero-normalized cross-correlation (ZNCC), the zero-normalized sum of squared differences (ZNSSD), or the parametric sum of squared differences (PSSD), is first defined to evaluate the similarity between the reference subset and its target subset. Then, proper image matching algorithms, such as the classic Newton– Raphson (NR) algorithm [6–9] and the state-of-the-art IC–GN algorithm [10], are then implemented to register the same physical points in these two images to retrieve the displacements and displacement gradients (strains).
∗
Corresponding author. E-mail address: [email protected] (L. Dong).
While displacement and strain fields can be directly obtained by using some of these DIC algorithms, the accuracy in strain results is proven to be low. Various approaches [11–23] were developed to obtain better strain calculations, such as the thin plate spline smoothing technique [11], the finite element method [12–14], the penalized least square regression method [15], the principal component analysis method [16], the fast Hermite element method [17], the local Kriging regression method [18], and the pointwise least squares (PLS) algorithm [19–21], etc. These methods can reduce the noise and the error of strains directly obtained by the NR or IC–GN algorithms, based on the idea of smoothing the displacement field in various ways, such as using RBF, Kriging, PCA, local polynomial, etc. Among the above-mentioned approaches, the PLS algorithm becomes widely accepted due to its effectiveness and simplicity. However, the PLS algorithm still has some limitations [24]. First, for inhomogeneous deformations, it is difficult to find a proper strain calculation window that can keep a balance between the effect of displacement-smoothing and the accuracy of strain-field calculation. Second, the accuracy and the robustness of the PLS algorithm is sensitive to the order of the polynomial kernel function, therefore it may give poor results for large-gradient or highly-oscillating strain fields when the order of polynomial does not correspond to the real deformation [25]. Besides, there are few studies about strain measurement of largegradient and highly-oscillating strain fields [26,27]. In the recent paper [27], the local Hermite method (LH) with C2 elements [27] proposed by Xin Li et al. is reported to generally produce lower error than the
https://doi.org/10.1016/j.optlaseng.2019.05.012 Received 13 July 2018; Received in revised form 8 May 2019; Accepted 8 May 2019 Available online 17 May 2019 0143-8166/© 2019 Elsevier Ltd. All rights reserved.
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
PLS algorithm in computing high-gradient strain fields. However, the LH method is roughly 10 times slower than the latter with the same window size [27]. In this work, we propose an alternative weak-form framework [28,29] to calculate strain fields from noisy displacement measurements by DIC. The regularly used subpixel registration algorithms are used to calculate displacements. Direct differentiation or smoothing the displacement field, however, is not used to calculate displacementgradients/strains. Inspired by the idea of the weak-form, which is frequently used in modern computational mechanics methods such as finite elements and boundary elements, the currently proposed method transfers differentiation into integration by applying the theorem of integration by parts. In order to derive a closed-form formula of strains, the cosine series is used as a trial as well as test functions to explore their orthogonality. Moreover, Filon’s method [30] is used to improve the accuracy of numerical integration as sinusoidal terms are involved. The effectiveness of the proposed method is verified by using numerical simulations. The results reveal that the proposed method is more accurate and less sensitive to noise when large-gradient or highly-oscillating strain fields are to be measured, and has similar computational efficiency as compared to the PLS algorithm in most of these numerical simulations. The rest of this paper is organized as follows. In Section 2, the IC–GN algorithm and the PLS algorithm for displacement and strain calculations are briefly reviewed. In Section 3, the derivation of the proposed weak-form method is given in detail. In Section 4, numerical experiments with various displacement/strain fields are used to demonstrate the advantages of the proposed method. In Section 5, we complete this paper with some concluding remarks.
where p denotes the unknown parameter vector containing displacements and displacement gradients of the subset. f(xi , yj ) and 𝑔(𝑥′𝑖 , 𝑦′𝑗 ) are the grayscale levels at coordinates (xi , yj ) in the reference subset and (𝑥′𝑖 , 𝑦′𝑗 ) in the target subset; 𝑓𝑚 =
𝑀 𝑀 ∑ ∑ 1 (2𝑀+1)2 𝑖=−𝑀 𝑖=−𝑀
𝑓 (𝑥𝑖 , 𝑦𝑗 ), and
𝑀 𝑀 ∑ ∑ 1 (2𝑀+1)2 𝑖=−𝑀 𝑖=−𝑀
𝑔(𝑥′ 𝑖 , 𝑦′ 𝑗 ) are the mean intensity values of the
𝑢(𝑥, 𝑦) = 𝑎0 + 𝑎1 𝑥 + 𝑎2 𝑦 𝑣(𝑥, 𝑦) = 𝑏0 + 𝑏1 𝑥 + 𝑏2 𝑦
(2)
𝑔𝑚 =
two subsets. The coordinate (𝑥′𝑖 , 𝑦′𝑗 ) is mapped from (xi , yj ) by displacement mapping function. In order to calculate the displacements and displacement gradients of the points of interest, the IC–GN algorithm is utilized to optimize the correlation criterion given in Eq. (1) in our work. The PLS algorithm for strain estimation is given as follows. A square strain calculation window containing (2𝑚 + 1) × (2𝑚 + 1) discrete points centered at the point of interest is chosen firstly. In the strain calculation window, the displacement field can be approximated as polynomial [18]. If a linear approximation is used, we have:
where x and y are the coordinates of discrete points within the strain calculation window; ai ,bi (𝑖 = 0, 1, 2) are the unknown polynomial coefficients, which can be determined by the least-squares method. The strains at the center point can be obtained based on the polynomial coefficients as it is a linear displacement field. 3. Strain field measurement using the weak-form method As described in [24], different strain calculation windows of PLS will lead to different results for an inhomogeneous deformation. When the calculation window is too small, the noise cannot be effectively removed. A large window on the other hand leads to over-smoothed strains. Moreover, errors will occur when a small calculation window is used to deal with large-gradient or highly oscillating strain fields, even if little or no noise is contained in displacement. The proposed weak-form method, which will overcome the shortcomings mentioned above, is described in detail in the following two subsections. In this paper, we only give the procedure to calculate the x-direction displacement-gradient/strain, while other strains can be calculated using similar algorithms.
2. Displacement and strain fields measurement using DIC DIC deals with the reference image and the deformed image of the test sample’s surface. The basic principle of the standard subset-based DIC is schematically illustrated in Fig. 1. A square reference subset of (2𝑀 + 1) × (2𝑀 + 1) pixels centered at the interrogated point P(x, y) is usually chosen to track its most similar target subset in the deformed image using certain correlation criteria. The Zero-mean Normalized Sum of Squared Differences (ZNSSD) criterion [6], which is insensitive to the scale and offset changes in illumination lighting fluctuations, is employed in this work:
⎡ ⎤ ( ) ( ) 𝑀 𝑀 ⎢ ⎥ ∑ ∑ 𝑓 𝑥𝑖 , 𝑦𝑗 − 𝑓𝑚 𝑔 𝑥′ 𝑖 , 𝑦′ 𝑗 − 𝑔𝑚 ⎢√ ⎥ 𝐶ZNSSD (𝑝) = −√ ⎢ ∑ [ ( ) ]2 [ ( ) ]2 ⎥ ∑𝑀 ∑𝑀 ∑𝑀 𝑖=−𝑀 𝑖=−𝑀 𝑀 ′ ′ ⎢ ⎥ 𝑖=−𝑀 𝑗=−𝑀 𝑓 𝑥𝑖 , 𝑦𝑗 − 𝑓𝑚 𝑖=−𝑀 𝑗=−𝑀 𝑔 𝑥 𝑖 , 𝑦 𝑗 − 𝑔𝑚 ⎦ ⎣
(1)
Fig. 1. The reference image and the deformed image used in the DIC method.
496
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
3.1. The weak-form method
3.2. Filon’s integration method
Consider the infinitesimal strain component in the x-direction: 𝜕𝑈 (𝑥, 𝑦) 𝜀𝑥 (𝑥, 𝑦) = , 0 ≤ 𝑥 ≤ 𝑎, 0 ≤ 𝑦 ≤ 𝑏, 𝜕𝑥
Consider an integral of this type: ( ) 𝑏 𝑘𝜋𝑦 𝐼 = ∫ cos 𝑈 (𝑎, 𝑦)𝑑𝑦 𝑏 0
(3)
where a and b are the length and the width of the rectangle region of interest, respectively; U(x, y) is the x-direction displacement obtained by DIC at discrete points. By multiplying both sides of Eq. (3) by a test function P(x, y) and applying the domain integration, we can obtain a weak-form integral equation: 𝑏 𝑎
𝑏 𝑎
∫ ∫ 𝜀𝑥 (𝑥, 𝑦)𝑃 (𝑥, 𝑦)𝑑 𝑥𝑑 𝑦 = ∫ ∫ 0 0
0 0
𝜕𝑈 (𝑥, 𝑦) 𝑃 (𝑥, 𝑦)𝑑 𝑥𝑑 𝑦 𝜕𝑥
Dividing the range (0, b) into 2n equal intervals. The interval length is: ℎ=𝑏∕2𝑛 And we denote: 𝑘𝜋 ℎ 𝑏
(4)
=𝜃 𝑟ℎ = 𝑦𝑟 𝑈 (𝑎, 𝑟ℎ) = 𝑈𝑟
(5)
where r is an integer. After a series of derivations, we have: ( ) [ ( ) ] 𝑏 𝑘𝜋𝑦 𝑘𝜋𝑏 + 𝛽𝑆2𝑟 + 𝛾𝑆2𝑟−1 𝑈 (𝑎, 𝑦)𝑑𝑦 = ℎ 𝛼𝑈 (𝑎, 𝑏) sin ∫ cos 𝑏 𝑏 0
Integrating the right side of Eq. (4) by parts we have: ∫0𝑏 ∫0𝑎 𝜀𝑥 (𝑥, 𝑦)𝑃 (𝑥, 𝑦)𝑑 𝑥𝑑 𝑦 = (𝑥,𝑦) 𝑑𝑥𝑑𝑦 ∫0𝑏 [𝑃 (𝑎, 𝑦)𝑈 (𝑎, 𝑦) − 𝑃 (0, 𝑦)𝑈 (0, 𝑦)]𝑑𝑦 − ∫0𝑏 ∫0𝑎 𝑈 (𝑥, 𝑦) 𝜕𝑃𝜕𝑥
It should be noted that, through integration by parts, the necessity of differentiating U is avoided, in the price of differentiating the test function P. In this way, we only need discrete values U when evaluating Eq. (4), so that the problem of numerical differentiation is transferred to numerical integration. In this work, ɛx (x, y) is expanded in Fourier series. Through expanding the region of interest to {(𝑥, 𝑦)|−𝑎 ≤ 𝑥 ≤ 𝑎, −𝑏 ≤ 𝑦 ≤ 𝑏} and imposing symmetry of ɛx (x, y) (i.e.𝜀𝑥 (𝑥, 𝑦) = 𝜀𝑥 (−𝑥, 𝑦) = 𝜀𝑥 (𝑥, −𝑦) = 𝜀𝑥 (−𝑥, −𝑦)), we can simplify the complete Fourier series as the cosine series: 𝜀𝑥 (𝑥, 𝑦) =
𝑁 ∑ 𝑀 ∑ 𝑛=0 𝑚=0
ℎ𝑛𝑚 cos
(
) ( 𝑚𝜋𝑦 ) 𝑛𝜋𝑥 cos 𝑎 𝑏
𝑆2𝑟 =
𝑟=0
( ) ( ) 𝑘𝜋𝑦2𝑟 𝑈 𝑎, 𝑦2𝑟 cos − 12 𝑈 (𝑎, 0) − 12 𝑈 (𝑎, 𝑏) cos 𝑘𝜋 𝑏
𝑛 ∑
𝑟=1
( ) 𝑈 𝑎, 𝑦2𝑟−1 cos
𝑘𝜋𝑦2𝑟−1 𝑏
(14)
when k is not equal to 0: 𝛼= 1∕𝜃+cos𝜃sin𝜃∕𝜃 2 − 2sin2 𝜃∕𝜃 3 𝛽= 2[(1 + cos2 𝜃)∕𝜃 2 − 2sin𝜃cos𝜃∕𝜃 3 ( ) 𝛾= 4 sin𝜃∕𝜃 3 − cos𝜃∕𝜃 2
(6)
when k is equal to 0: 𝛼= 0 𝛽= 2∕3 𝛾= 4∕3 Details about Filon’s integration can be found in [30]. 4. Verification using numerical experiments To investigate the performance of the weak-form method proposed in this paper, five cases of numerical experiments were conducted in this study: (1) a small homogeneous deformation, (2) two sinusoidal displacement fields, (3) a Gaussian displacement field deformations, (4) an analogous sinusoidal-Gaussian displacement field and (5) deformation in a perforated plate. In the case (1–4), computer-simulated images are used and the algorithm of generating computer-simulated speckle images is given in [31,32]. The simulated images adopt the following parameters: the size of simulated speckle images is 601 × 601 pixels, the number of speckle granules is 13,000 and the size of speckles is 2.5 pixels. Gaussian noises are added to both reference and deformed images to simulate the image noises encountered in practical experiments. The displacements, which are at regularly distributed 17,161 (131 × 131) points with a grid step of 4 pixels for the case (1) and 75,625 (275 × 275) points with a grid step of 2 pixels for the case (2–4), are computed using the IC–GN algorithm implemented in the software Ncorr taken from the website www.ncorr.com. In the case (5), the noisy displacement is created by adding normally distributed random numbers to the analytical displacement field. Then the proposed weak-form method is used to calculate the strain fields. In order to compare the accuracy of the proposed weak-form method with the LH method, we also use the LH method with C2 elements to calculate the strain fields of cases (3–5). All parameters of the LH method adopted in this paper can be found in [27]. To quantitatively investigate the performance of the proposed weakform method, the root-mean-square (RMS) error is calculated. The RMS error is defined as follows:
(8)
Substituting Eq. (6) into Eq. (8), and by using the orthogonality: ( ) ( ) 𝑗𝜋𝑦 𝑘𝜋𝑦 𝑏 cos 𝑑𝑦 = 𝛿𝑗𝑘 (9) ∫ cos 𝑏 𝑏 2 0
[ ( ) ( ) ] 𝑒𝑛𝑚 = ∫0𝑏 cos (𝑛𝜋) cos 𝑚𝜋𝑦 𝑈 (𝑎, 𝑦) − cos 𝑚𝜋𝑦 𝑈 (0, 𝑦) 𝑑𝑦+ 𝑏) 𝑏 ( ( ) sin 𝑛𝜋𝑥 cos 𝑚𝜋𝑦 𝑑 𝑥𝑑 𝑦 ∫0𝑏 ∫0𝑎 𝑈 (𝑥, 𝑦) 𝑛𝜋 𝑎 𝑎 𝑏
𝑛 ∑
𝑆2𝑟−1 =
𝑏
we have: 4𝜆 𝑛 ℎ𝑛𝑚 = 𝑒 𝑎𝑏 𝑚
(13)
where
where ℎ𝑛𝑚 are unknown coefficients. And (𝑁 + 1) × (𝑀 + 1) is called the number of expanding terms of the cosine series in this paper. According to the numerical results and our experience, we suggest that the range 𝐼 𝐽 of expanding terms is 15 × 15 ≤ (𝑁 + 1) × (𝑀 + 1) ≤ 𝐼3 × 𝐽3 , where I × J is the number of points calculated in the region of interest. To get a closedform expression of ℎ𝑛𝑚 , we utilize the orthogonality of cosine functions and choose test function P(x, y) as follows: ( ) ( ) 𝑗𝜋𝑥 𝑘𝜋𝑦 𝑃 (𝑥, 𝑦) = cos cos (7) 𝑎 𝑏 Substituting P(x, y) into Eq. (5), we have: ( ) ( ) cos 𝑘𝜋𝑦 𝑑 𝑥𝑑 𝑦 ∫0𝑏 ∫0𝑎 𝜀𝑥 (𝑥, 𝑦) cos 𝑗𝜋𝑥 𝑎 𝑏 [ ( ) ( ) ] 𝑘𝜋𝑦 𝑏 = ∫0 cos (𝑗𝜋) cos 𝑏 𝑈 (𝑎, 𝑦) − cos 𝑘𝜋𝑦 𝑈 (0, 𝑦) 𝑑𝑦+ 𝑏 ( ) ( ) cos 𝑘𝜋𝑦 𝑑 𝑥𝑑 𝑦 sin 𝑗𝜋𝑥 ∫0𝑏 ∫0𝑎 𝑈 (𝑥, 𝑦) 𝑗𝜋 𝑎 𝑎 𝑏
(12)
(10)
(11)
where, 𝜆 = 1∕4 when 𝑚 = 0 and 𝑛 = 0; 𝜆 = 1∕2 when one of m and n is 0; 𝜆 = 1 when both of m and n is not equal to 0. Strains distributions can be calculated by directly evaluating the integral of Eq. (11). However, traditional numerical integration algorithms, like the trapezoidal rule and the Gaussian integration, will give poor results, due to the sinusoidal terms in the integrand. Thus, Filon’s integration method is used which are specifically designed for integrands with an oscillatory factor of sin (kx)orcos (kx). 497
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 2. The x-directional displacement filed obtained by the IC–GN algorithm for the case of homogeneous deformation.
√ 𝑅𝑀𝑆 =
∑
(𝜀𝑚 − 𝜀)2 𝑛
(15)
where ɛm denotes the strain calculated by the PLS algorithm, or the proposed weak-form method at each calculation point; ɛ denotes the exact value of strains; n is the number of calculation points. 4.1. Homogeneous deformation As the simplest case of deformations with small strain-gradients, a homogeneous deformation of uniaxial tension is applied in x-direction: 𝑈 = 0.002𝑥, 𝜀𝑥 = 0.002, 0 ≤ 𝑥 ≤ 600, while the displacement V in y-direction is set to 0. Gaussian noise with a mean value of zero and a standard deviation of 0.5% of the entire 256 grayscales is added to both reference and deformed images. Fig. 2 shows the displacement field U obtained by the IC–GN algorithm. To obtain a reliable strain field, the proposed weak-form method is used to post-process the displacement calculated by the IC–GN algorithm. Fig. 3(a) gives the strain field obtained by the weak-form method with 9 × 9expanding terms (the RMS error is 98.36𝜇ɛ). For comparison, the strain field calculated using the PLS algorithm with the best strain calculation window of 21 × 21points (the RMS error is 70.93𝜇ɛ) is also shown in Fig. 3(b). Further, the effects of expanding terms of the proposed weakform method and the calculation windows size of the PLS algorithm are investigated. Fig. 4 shows the RMS errors of strain fields calculated by the proposed weak-form method with expanding terms ranging from 9 × 9 to 18 × 18, and by the PLS algorithm with various strain calculation windows varying from 3 × 3 points to 21 × 21 points. It can be found that the error dispersion of the weak-form method with 10 different expanding terms is smaller than that of the PLS algorithm with 10 different strain calculation windows. For the homogeneous deformation, the small number of expanding terms for the weak-form method will give a more accurate result, while large strain calculation windows are good choices for the PLS algorithm.
Fig. 3. The x-directional strain field obtained by (a) the weak-form method, and (b) the PLS algorithm, for the case of homogeneous deformation.
Example 1: The imposed displacement field of the first example is as follows: 𝑈 = sin
(
) ( ) 2𝜋 2𝜋𝑥 2𝜋𝑥 , 𝜀𝑥 = , 0 ≤ 𝑥 ≤ 600, cos 300 300 300
Gaussian noise with a mean value of zero and a standard deviation of 1% of the entire 256 grayscales is also added to both reference and deformed images. Fig. 5 shows the displacement field obtained by the IC–GN algorithm. Fig. 6 shows the strain field calculated by different algorithms. The side view of the strain field obtained by using the weak-form method with18 × 18 expanding terms (the RMS error is629.69𝜇ɛ) is shown in Fig. 6(a). Also, the side view of the strain field calculated using the PLS algorithm with the best strain window of 11 × 11 points (the RMS error is758.26𝜇ɛ) is shown in Fig. 6(b). For comparison, the side view of the exact strain field is shown in Fig. 6(c). Further, the effects of expanding terms of the proposed weak-form method and the strain calculation windows size of the PLS algorithm are investigated. Fig. 7 shows the RMS errors of strain fields calculated by the proposed weak-form method with expanding terms ranging from 18 × 18 to 27 × 27, and by the PLS algorithm with various strain calculation windows varying from 3 × 3 points to 21 × 21 points.
4.2. Sinusoidal displacement fields To study the performance of the proposed weak-form method for large strain-gradient or highly-oscillatory deformations, the next two examples are considering in-plane sinusoidal displacement fields in the x-direction. 498
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 4. RMS errors of the estimated strain fields with various expanding terms of the weak-form method and various strain calculation windows of the PLS algorithm, for the case of homogeneous deformation.
Fig. 5. The x-directional displacement field obtained by the IC–GN algorithm for example 1 of the sinusoidal displacement fields.
Example 2: The second example with larger strain gradients and larger noise than the first example is given as follows: ( ) ( ) 2𝜋 2𝜋𝑥 2𝜋𝑥 𝑈 = sin cos , 𝜀𝑥 = , 0 ≤ 𝑥 ≤ 600, 200 200 200 Gaussian noise with a mean value of zero and a standard deviation of 2% of the entire 256 grayscales is also added to both reference and deformed images. Fig. 8 shows the computed x-direction displacement field using the IC–GN algorithm. By post-processing the displacement calculated by the IC–GN algorithm using the weak-form method with 18 × 18 expanding terms, the side view of the strain field is shown in Fig. 9(a). Besides, the side view of the strain field calculated using the PLS algorithm with the best strain window of 11 × 11 points is shown in Fig. 9(b). For comparison, the side view of the exact strain field is shown in Fig. 9(c). In a similar manner, the errors of strain fields corresponding to various expanding terms of the weak-form method and different strain calculation windows of the PLS algorithm are shown in Fig. 10. From these two examples, it is seen that the proposed method has advantages for deformations with large strain-gradients. The first ad-
Fig. 6. Side views of strain fields obtained by (a) the weak-form method, (b) the PLS algorithm, and (c) the exact strain, for example 1 of the sinusoidal displacement fields. 499
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 7. RMS errors of the estimated strain fields with various expanding terms of the weak-form method and various strain calculation windows of the PLS algorithm, for example 1 of the sinusoidal displacement fields.
Fig. 8. The x-directional displacement field obtained by the IC–GN algorithm for example 2 of the sinusoidal displacement fields.
vantage of the weak-form method is that the error dispersion of strains with 10 different expanding terms is smaller than the PLS algorithm with 10 different strain calculation windows. The other advantage is that strain fields calculated by the weak-form method with various expanding terms are more accurate and less insensitive to noise. 4.3. A Gaussian displacement field A Gaussian displacement field proposed by Xu et al. [33] is considered here: ( ) √ 𝑈 = 𝑎 × 𝑐 𝜋∕2𝑒𝑟𝑓 √𝑥 , −300 ≤ 𝑥 ≤ 300 2𝑐
𝜀𝑥 = 𝑎
2 −𝑥 × 𝑒 2𝑐 2
where parameter a is the intensity of the Gauss function and c is the standard error of the Gauss function. In this paper, a is 0.0025 and c is 16. Gaussian noise with a mean value of zero and a standard deviation of 1% of the entire 256 grayscales is added to both reference and deformed images.
Fig. 9. Side views of strain fields obtained by (a) the weak-form method, (b) the PLS algorithm, (c) the exact strain, for example 2 of the sinusoidal displacement fields. 500
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 10. RMS errors of the estimated strain fields with various expanding terms of the weak-form method and various strain calculation windows of the PLS algorithm, for example 2 of the sinusoidal displacement fields.
Fig. 11. The side view of the x-directional displacement field obtained by the IC–GN algorithm for the Gaussian displacement field.
Fig. 11 shows the computed x-direction displacement field using the IC–GN algorithm. By post-processing the obtained displacement using the weak-form method with 36 × 36 expanding terms, the side view of the strain field is shown in Fig. 12(a). Besides, the side views of the strain field calculated using the PLS algorithm with the best strain window of 7 × 7 points and using the LH method with the best strain window of 29 × 29 points are shown in Fig. 12(b) and Fig. 12(c), respectively. The errors of strain fields corresponding to various expanding terms of the weak-form method and different strain calculation windows of the PLS algorithm as well as the LH method are shown in Fig. 13. 4.4. An analogous sinusoidal-Gaussian displacement field An analogous sinusoidal-Gaussian displacement field proposed by Yang et al. [34] is considered here: 𝑈 (𝑥, 𝑦) = sin (2𝜋𝐸 (𝑥)) sin (2𝜋𝐸 (𝑦))∕10, ) 2 ( 𝐸 (𝑥) = 𝑒−(𝑥−300) ∕ 2𝜎 2 ) 2 ( − ( 𝑦 −300) 2 𝐸 (𝑦 ) = 𝑒 ∕ 2𝜎
0 ≤ 𝑥 ≤ 600 Fig. 12. Side views of strain fields obtained by (a) the weak-form method, (b) the PLS algorithm, (c) the LH method, and (d) the exact strain, for the Gaussian displacement field. 501
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 13. RMS errors of the estimated strain fields with various expanding terms of the weak-form method, various strain calculation windows of the PLS algorithm and the LH method, for the Gaussian displacement field.
where 𝜎 denotes the Gaussian Root-Mean-Square (RMS) width, which is 80 in this paper. Gaussian noise with a mean value of zero and a standard deviation of 1% is added to both reference and deformed images. Fig. 14 shows the side views of the computed x-direction displacement field using the IC–GN algorithm. By post-processing the displacement using the weak-form method with 27 × 27 expanding terms, the side view of the strain field is shown in Fig. 15(a). Besides, the side views of strain field calculated using the PLS algorithm with the best strain window of 7 × 7 points and using the LH method with the best strain window of 29 × 29 points are shown in Fig. 15(b) and Fig. 15(c), respectively. The errors of strain fields corresponding to a various number of expanding terms of the weak-form method and different strain calculation windows of the PLS algorithm as well as the LH method are shown in Fig. 16.
4.5. Deformation in a perforated plate To validate the feasibility of the proposed method for perforated plate problems, the last case investigated is a plate with a small hole under a uniform tension in x-direction, as shown in Fig. 17. The material properties of the plate are 𝐸 = 0.8𝐺𝑃 𝑎 and 𝑣 = 0.3. The radius of the hole is 50mm and the tension q is 23 × 107 𝑁∕𝑚2 . The size of the region of interest is 598 × 598mm2 with a grid step of 2mm. The analytical solutions of the displacement field and the strain field can be found in [35], and the noisy displacement is created by adding normally distributed random numbers to the analytical displacement field. In this example, the weak-form method described in Section 3.1 cannot be directly used, because of the existence of the hole in the center of the plate. A simple pre-processing of the displacement data is implemented before calculating the strain field: dividing the original image into 4 subdomains, and extrapolating displacements in each subdomain to a rectangular box. Fig. 18(a) illustrates an irregular region of interest, in which displacements at red points are already calculated by NR, IC–GN or other algorithms. In Fig. 18(b), the virtual displacements at blue points are to be extrapolated by using displacements at red points. The dataextrapolating process is described as follows: (1) a square calculation window containing 5 × 5 − 1 red points and 1 blue point (point P) is chosen to calculate the virtual displacement at the point P, as shown in Fig. 18(b). (2) The displacement field in the calculation window can
Fig. 14. (a)Theoretical x-directional displacement in 3D view, (b) theoretical x-directional displacement in side view, and (c) x-directional displacements obtained by the IC–GN algorithm in side view.
502
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 16. RMS errors of the estimated strain fields with various expanding terms of the weak-form method, various strain calculation windows of the PLS algorithm and the LH method, for the analogous sinusoidal-Gaussian displacement field.
Fig. 17. A plate with a small hole submitted to a uniform tension in x-direction.
be approximated by using the polynomial 𝑈 = 𝑎0 + 𝑎1 𝑥 + 𝑎2 𝑦 + 𝑎3 𝑥2 + 𝑎4 𝑥𝑦 + 𝑎5 𝑦2 , in which the unknown polynomial coefficientsa0 , a1 , a2 , a3 , a4 , a5 are solved by the least-squares method using displacements at blue points. (3) After the virtual displacement at the point P is obtained, it will be treated as a known value to calculate the virtual displacement at the next blue point P1 . The same procedure is used to calculate virtual displacements at each blue point from left to right and from top to bottom in turn. Finally, when all virtual displacements at blue points are calculated, the proposed weak-form method can be used to calculate the strain field in the rectangular region shown in Fig. 18(b). For the perforated plate problem, it should be noted that the order of polynomial using in extrapolating process will affect the accuracy of the strain to some extent. By dividing the region of interest into 4 subdomains as shown in Fig. 17, the strain field can be calculated by using the process described above in each subdomain. Fig. 19 shows ɛx calculated using the proposed weak-form method with 30 × 30expanding terms in each subdomain. Fig. 20 shows the ɛx along the y-axis direction calculated by the proposed weak-form method, the PLS algorithm, and the analytical solution, respectively. Fig. 21 displays the RMS error of strain field corresponding to various expanding terms of the weak-form method, different strain calculation windows of the PLS algorithm and the LH method. It is observed that the weak-form method can improve the accuracy of
Fig. 15. Side views of strain fields obtained by (a) the weak-form method, (b) the PLS algorithm, and (c) the LH method for the analogous sinusoidal-Gaussian displacement field.
503
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 19. The x-directional strain field of the perorated plate submitted to a uniform tension obtained by the weak-form method.
Fig. 20. The x-directional strain field along the y-axis of the perorated plate obtained by the weak-form method, the PLS algorithm with the strain window of 5 × 5 as well as 11 × 11 points, and the exact value.
Fig. 18. (a) an irregular region of interest, (b) a rectangle region by extrapolating displacements for the irregular region of interest. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)
strain fields especially in obvious stress concentrated areas compared to the PLS algorithm. In order to elucidate the computational efficiency of the proposed method, the computational time of above all numerical experiments on the condition of the best-expanding terms and the best strain window is shown in Table 1. It is seen that the proposed method has similar computational efficiency as compared to the PLS algorithm in most of these numerical simulations.
Fig. 21. RMS errors of the estimated strain fields with various expanding terms of the weak-form method, various strain calculation windows of the PLS algorithm and the LH method, for the perorated plate. 504
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Table 1 The computational time for all numerical experiments on the condition of the best expanding terms of the weak-form method and the best strain windows of the PLS algorithm. EF is the expanding terms of the weak-form method. SCW is the strain calculation windows of the PLS algorithm. Method time (ET or SCW) Numerical experiments
The weak-form method
The PLS algorithm
Homogeneous deformation Example 1 of the sinusoidal displacement field Example 2 of the sinusoidal displacement field The Gaussian displacement field The analogous sinusoidal-Gaussian displacement field The deformation in a perforated plate
0.33 s (ET:9 × 9) 4.02 s (ET:18 × 18) 4.23 s (ET:18 × 18) 10.71 s (ET:29 × 29) 9.10 s (ET:27 × 27) 11.04 s (ET:30 × 30)
14.31 s (SCW:21 × 21) 20.29 s (SCW:11 × 11) 20.71 s (SCW:11 × 11) 11.34 s (SCW:7 × 7) 11.23 s (SCW:7 × 7) 8.00 s (SCW:5 × 5)
5. Conclusion
[7] Bruck HA, McNeill SR, Sutton MA, Peters WH. Digital image correlation using Newton–Raphson method of partial differential correction. Exp Mech 1989;29:261–7. [8] Lu H, Cary PD. Deformation measurements by digital image correlation: implementation of a second-order displacement gradient. Exp Mech 2000;40:393–400. [9] Yoneyama S, Morimoto Y. Accurate displacement measurement by correlation of colored random patterns. JSME Int J Series A 2003;46:178–84. [10] Baker S, Matthews I. Lucas–Kanade 20 years on: a unifying framework. Int J Comput Vis 2004;56:221–55. [11] Wang CCB, Deng JM, Ateshian GA, Hung CT. An automated approach for direct measurement of two-dimensional strain distributions within articular cartilage under unconfined compression. J Biomech Eng 2002;124:557–67. [12] Sutton MA, Turner JL, Bruck HA, Chae TA. Full-field representation of discretely sampled surface deformation for displacement and strain analysis. Exp Mech 1991;31:168–77. [13] Tong W. Detection of plastic deformation patterns in a binary aluminum alloy. Exp Mech 1997;37:452–9. [14] Yoneyama S. Smoothing measured displacements and computing strains utilising finite element method. Strain 2011;47:258–66. [15] Han Y, Kim DW, Kwon HJ. Application of digital image cross-correlation and smoothing function to the diagnosis of breast cancer. J Mech Behav Biomed Mater 2012;14:7–18. [16] Grama SN, Subramanian SJ. Computation of full-field strains using principal component analysis. Exp Mech 2014;54:913–33. [17] Zhao J, Song Y, Wu X. Fast Hermite element method for smoothing and differentiating noisy displacement field in digital image correlation. Opt Lasers Eng 2015;68:25–34. [18] Wang D, DiazDelaO FA, Wang W, Lin X, Patterson EA, Mottershead JE. Uncertainty quantification in DIC with Kriging regression. Opt Lasers Eng 2016;78:182–95. [19] Wattrisse B, Chrysochoos A, Muracciole J-M, Némoz-Gaillard M. Analysis of strain localization during tensile tests by digital image correlation. Exp Mech 2001;41:29–39. [20] Pan B, Xie H, Guo Z, Hua T. Full-field strain measurement using a two-dimensional Savitzky–Golay digital differentiator in digital image correlation. Opt Eng 2007;46:033601. [21] Pan B, Asundi A, Xie H, Gao J. Digital image correlation using iterative least squares and pointwise least squares for displacement field and strain field measurements. Opt Lasers Eng 2009;47:865–74. [22] Lehoucq RB, Reu PL, Turner DZ. A novel class of strain measures for digital image correlation. Strain 2015;51:265–75. [23] Dai X, Yang F, Chen Z, Shao X, He X. Strain field estimation based on digital image correlation and radial basis function. Opt Lasers Eng 2015;65:64–72. [24] Pan B, Yuan J, Xia Y. Strain field denoising for digital image correlation using a regularized cost-function. Opt Lasers Eng 2015;65:9–17. [25] Schreier HW, Sutton MA. Systematic errors in digital image correlation due to undermatched subset shape functions. Exp Mech 2002;42:303–10. [26] Hwang SF, Wu WJ. Deformation measurement around a high strain-gradient region using a digital image correlation method. J Mech Sci Technol 2012;26:3169–75. [27] Li X, Fang G, Zhao JQ, Zhang ZM, Wu XX. Local Hermite (LH) Method: an accurate and robust smooth technique for high-gradient strain reconstruction in digital image correlation. Opt Lasers Eng 2019;112:26–38. [28] Liu CS, Dong L. Closed-form higher-order numerical differentiators for differentiating noisy signals. Appl Math Comput 2019 in press. [29] Dong L, Alotaibi A, Mohiuddine SA, Atluri SN. Computational methods in engineering: a variety of primal & mixed methods, with global & local interpolations, for well-posed or ill-posed BCs. Comput Model Eng Sci 2014;99:1–85. [30] Filon LNG III. On a quadrature formula for trigonometric integrals. In: Proceedings of the Royal Society of Edinburgh, 49; 1930. p. 38–47. [31] Zhou P, Goodson KE. Subpixel displacement and deformation gradient measurement using digital image/speckle correlation (DISC). Opt Eng 2001;40:1613–20. [32] Pan B, Xie HM, Xu BQ, Dai FL. Performance of sub-pixel registration algorithms in digital image correlation. Measur Sci Technol 2006;17:1615–21. [33] Xu X, Su Y, Cai Y, Cheng T, Zhang Q. Effects of various shape functions and subset size in local deformation measurements using DIC. Exp Mech 2015;55:1575–90. [34] Yang X, Chen X, Xi J. Comparative analysis of warp function for digital image correlation-based accurate single-shot 3D shape measurement. Sensors 2018;18:1208. [35] Timoshenko SP, Goodier JN. Theory of elasticity; 1970.
The assessment of strain fields is important for deformation measurement. However, the strain fields obtained directly by the NR algorithm, the IC–GN algorithm, or genetic algorithms, have poor accuracy. The PLS algorithm is widely used for retrieving strain distribution. However, the PLS algorithm has some limitations when dealing with large gradients strain fields, especially with noisy measurements. In this paper, an alternative strain estimation approach based on the weak form is proposed, and verified by numerical simulations. Numerical simulations reveal that when strain gradients are small, strain fields obtained by the weak-form method are comparable to the PLS algorithm. However, for large strain-gradient or highly-oscillating deformations, the results obtained by the weak-form method are more accurate than the PLS algorithm. Furthermore, the error dispersion of the weak-form method is smaller than the PLS algorithm no matter strain gradients are small or large. For complex large strain-gradient or highlyoscillating deformation problems, although the LH method can obtain more accurate strain fields than the PLS algorithm, it is much more computationally expensive and complex than the weak-form method and the PLS algorithm. Through numerical experiments, one can also see that the number of expanding terms affects the accuracy of the weak-form method, although in a certain range the error dispersion is quite small. In the 𝐼 𝐽 suggested range [ 15 × 15 , 𝐼3 × 𝐽3 ] for the number of expanding terms, a small number is preferred for small strain-gradients and for large noises. It should also be noted that, while the weak-form method described in Section 3.1 is directly used to calculate the strain distribution for a rectangular region, irregular regions of interest can also be considered by extrapolating the original displacement data to a rectangular box. This is demonstrated in Section 4.5 with the numerical simulation of a perforated plate. Acknowledgment The authors thankfully acknowledge the support from the National Key Research and Development Program of China (No. 2017YFA0207800). References [1] Nakadate S, Yatagai T, Saito H. Electronic speckle pattern interferometry using digital image processing techniques. Appl Opt 1980;19:1879–83. [2] Bing P, Kemao Q, Huimin X, Anand A. Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review. Measur Sci Technol 2009;20:062001. [3] Pan B. Recent progress in digital image correlation. Exp Mech 2011;51:1223–35. [4] Peters WH, Ranson WF. Digital imaging techniques in experimental stress analysis. SPIE; 1982. p. 5. [5] Peters WH, Ranson WF, Sutton MA, Chu TC, Anderson J. Application of digital correlation methods to rigid body mechanics. SPIE; 1983. p. 5. [6] Pan B. Reliability-guided digital image correlation for image deformation measurement. Appl Opt 2009;48:1535–42.
505
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
A weak-form method for strain field calculation with digital image correlation Leiting Dong a,∗, Qiuhong Yang b, Junbo Wang a a b
School of Aeronautic Science and Engineering, Beihang University, Beijing 100191, China COMAC Shanghai Aircraft Design and Research Institute, Shanghai 201210, China
a r t i c l e
i n f o
Keywords: Digital image correlation Strain field measurement Denoising Weak-form Test functions Trial functions
a b s t r a c t In digital image correlation (DIC), iterative spatial-domain cross-correlation algorithms have been routinely used to extract displacement fields from the recorded images. Then, strain fields are computed from the noisy displacement fields by a proper numerical different approach, like pointwise least squares fitting, which is tricky in choosing optimal calculation parameters. This work proposes an alternative weak-form framework to calculate strain fields from noisy displacement measurements by DIC. While direct numerical differentiation may amplify the noise in the displacement fields, the proposed method transfers differentiation into integration by employing the idea of the weak-form, which is the foundation of modern computational mechanics. By selecting cosine series as trial/test functions and employing their orthogonality, a semi-analytical formula for calculating the strain field is given. Moreover, the Filon’s method is used in this study to improve the accuracy of numerical integration as sinusoidal terms are involved. The effectiveness of the proposed method is verified by numerical examples. It is found that the proposed method is comparable to the famous PLS algorithm when it is for simple strain fields, while it is better in accuracy and its insensitivity to noise when large-gradient or highly-oscillating strain fields are to be measured.
1. Introduction Surface deformation measurement of materials and structures is an important task of experimental mechanics. Various optical methods, such as holography interferometry, electronic speckle pattern interferometry (ESPI) [1], digital image correlation (DIC) [2,3], have been proposed and advocated for non-contact and full-field deformation measurement. Among these methods, DIC techniques proposed by Peters et al. [4,5] have been extensively investigated and used due to its simplicity, robustness and wide range of applicability [2]. The basic principle of DIC method is simple and straightforward, namely matching the subsets of the surface in the reference image and the deformed image to calculate the displacement fields and strain fields. To this purpose, certain correlation criteria [2], e.g., the zero-normalized cross-correlation (ZNCC), the zero-normalized sum of squared differences (ZNSSD), or the parametric sum of squared differences (PSSD), is first defined to evaluate the similarity between the reference subset and its target subset. Then, proper image matching algorithms, such as the classic Newton– Raphson (NR) algorithm [6–9] and the state-of-the-art IC–GN algorithm [10], are then implemented to register the same physical points in these two images to retrieve the displacements and displacement gradients (strains).
∗
Corresponding author. E-mail address: [email protected] (L. Dong).
While displacement and strain fields can be directly obtained by using some of these DIC algorithms, the accuracy in strain results is proven to be low. Various approaches [11–23] were developed to obtain better strain calculations, such as the thin plate spline smoothing technique [11], the finite element method [12–14], the penalized least square regression method [15], the principal component analysis method [16], the fast Hermite element method [17], the local Kriging regression method [18], and the pointwise least squares (PLS) algorithm [19–21], etc. These methods can reduce the noise and the error of strains directly obtained by the NR or IC–GN algorithms, based on the idea of smoothing the displacement field in various ways, such as using RBF, Kriging, PCA, local polynomial, etc. Among the above-mentioned approaches, the PLS algorithm becomes widely accepted due to its effectiveness and simplicity. However, the PLS algorithm still has some limitations [24]. First, for inhomogeneous deformations, it is difficult to find a proper strain calculation window that can keep a balance between the effect of displacement-smoothing and the accuracy of strain-field calculation. Second, the accuracy and the robustness of the PLS algorithm is sensitive to the order of the polynomial kernel function, therefore it may give poor results for large-gradient or highly-oscillating strain fields when the order of polynomial does not correspond to the real deformation [25]. Besides, there are few studies about strain measurement of largegradient and highly-oscillating strain fields [26,27]. In the recent paper [27], the local Hermite method (LH) with C2 elements [27] proposed by Xin Li et al. is reported to generally produce lower error than the
https://doi.org/10.1016/j.optlaseng.2019.05.012 Received 13 July 2018; Received in revised form 8 May 2019; Accepted 8 May 2019 Available online 17 May 2019 0143-8166/© 2019 Elsevier Ltd. All rights reserved.
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
PLS algorithm in computing high-gradient strain fields. However, the LH method is roughly 10 times slower than the latter with the same window size [27]. In this work, we propose an alternative weak-form framework [28,29] to calculate strain fields from noisy displacement measurements by DIC. The regularly used subpixel registration algorithms are used to calculate displacements. Direct differentiation or smoothing the displacement field, however, is not used to calculate displacementgradients/strains. Inspired by the idea of the weak-form, which is frequently used in modern computational mechanics methods such as finite elements and boundary elements, the currently proposed method transfers differentiation into integration by applying the theorem of integration by parts. In order to derive a closed-form formula of strains, the cosine series is used as a trial as well as test functions to explore their orthogonality. Moreover, Filon’s method [30] is used to improve the accuracy of numerical integration as sinusoidal terms are involved. The effectiveness of the proposed method is verified by using numerical simulations. The results reveal that the proposed method is more accurate and less sensitive to noise when large-gradient or highly-oscillating strain fields are to be measured, and has similar computational efficiency as compared to the PLS algorithm in most of these numerical simulations. The rest of this paper is organized as follows. In Section 2, the IC–GN algorithm and the PLS algorithm for displacement and strain calculations are briefly reviewed. In Section 3, the derivation of the proposed weak-form method is given in detail. In Section 4, numerical experiments with various displacement/strain fields are used to demonstrate the advantages of the proposed method. In Section 5, we complete this paper with some concluding remarks.
where p denotes the unknown parameter vector containing displacements and displacement gradients of the subset. f(xi , yj ) and 𝑔(𝑥′𝑖 , 𝑦′𝑗 ) are the grayscale levels at coordinates (xi , yj ) in the reference subset and (𝑥′𝑖 , 𝑦′𝑗 ) in the target subset; 𝑓𝑚 =
𝑀 𝑀 ∑ ∑ 1 (2𝑀+1)2 𝑖=−𝑀 𝑖=−𝑀
𝑓 (𝑥𝑖 , 𝑦𝑗 ), and
𝑀 𝑀 ∑ ∑ 1 (2𝑀+1)2 𝑖=−𝑀 𝑖=−𝑀
𝑔(𝑥′ 𝑖 , 𝑦′ 𝑗 ) are the mean intensity values of the
𝑢(𝑥, 𝑦) = 𝑎0 + 𝑎1 𝑥 + 𝑎2 𝑦 𝑣(𝑥, 𝑦) = 𝑏0 + 𝑏1 𝑥 + 𝑏2 𝑦
(2)
𝑔𝑚 =
two subsets. The coordinate (𝑥′𝑖 , 𝑦′𝑗 ) is mapped from (xi , yj ) by displacement mapping function. In order to calculate the displacements and displacement gradients of the points of interest, the IC–GN algorithm is utilized to optimize the correlation criterion given in Eq. (1) in our work. The PLS algorithm for strain estimation is given as follows. A square strain calculation window containing (2𝑚 + 1) × (2𝑚 + 1) discrete points centered at the point of interest is chosen firstly. In the strain calculation window, the displacement field can be approximated as polynomial [18]. If a linear approximation is used, we have:
where x and y are the coordinates of discrete points within the strain calculation window; ai ,bi (𝑖 = 0, 1, 2) are the unknown polynomial coefficients, which can be determined by the least-squares method. The strains at the center point can be obtained based on the polynomial coefficients as it is a linear displacement field. 3. Strain field measurement using the weak-form method As described in [24], different strain calculation windows of PLS will lead to different results for an inhomogeneous deformation. When the calculation window is too small, the noise cannot be effectively removed. A large window on the other hand leads to over-smoothed strains. Moreover, errors will occur when a small calculation window is used to deal with large-gradient or highly oscillating strain fields, even if little or no noise is contained in displacement. The proposed weak-form method, which will overcome the shortcomings mentioned above, is described in detail in the following two subsections. In this paper, we only give the procedure to calculate the x-direction displacement-gradient/strain, while other strains can be calculated using similar algorithms.
2. Displacement and strain fields measurement using DIC DIC deals with the reference image and the deformed image of the test sample’s surface. The basic principle of the standard subset-based DIC is schematically illustrated in Fig. 1. A square reference subset of (2𝑀 + 1) × (2𝑀 + 1) pixels centered at the interrogated point P(x, y) is usually chosen to track its most similar target subset in the deformed image using certain correlation criteria. The Zero-mean Normalized Sum of Squared Differences (ZNSSD) criterion [6], which is insensitive to the scale and offset changes in illumination lighting fluctuations, is employed in this work:
⎡ ⎤ ( ) ( ) 𝑀 𝑀 ⎢ ⎥ ∑ ∑ 𝑓 𝑥𝑖 , 𝑦𝑗 − 𝑓𝑚 𝑔 𝑥′ 𝑖 , 𝑦′ 𝑗 − 𝑔𝑚 ⎢√ ⎥ 𝐶ZNSSD (𝑝) = −√ ⎢ ∑ [ ( ) ]2 [ ( ) ]2 ⎥ ∑𝑀 ∑𝑀 ∑𝑀 𝑖=−𝑀 𝑖=−𝑀 𝑀 ′ ′ ⎢ ⎥ 𝑖=−𝑀 𝑗=−𝑀 𝑓 𝑥𝑖 , 𝑦𝑗 − 𝑓𝑚 𝑖=−𝑀 𝑗=−𝑀 𝑔 𝑥 𝑖 , 𝑦 𝑗 − 𝑔𝑚 ⎦ ⎣
(1)
Fig. 1. The reference image and the deformed image used in the DIC method.
496
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
3.1. The weak-form method
3.2. Filon’s integration method
Consider the infinitesimal strain component in the x-direction: 𝜕𝑈 (𝑥, 𝑦) 𝜀𝑥 (𝑥, 𝑦) = , 0 ≤ 𝑥 ≤ 𝑎, 0 ≤ 𝑦 ≤ 𝑏, 𝜕𝑥
Consider an integral of this type: ( ) 𝑏 𝑘𝜋𝑦 𝐼 = ∫ cos 𝑈 (𝑎, 𝑦)𝑑𝑦 𝑏 0
(3)
where a and b are the length and the width of the rectangle region of interest, respectively; U(x, y) is the x-direction displacement obtained by DIC at discrete points. By multiplying both sides of Eq. (3) by a test function P(x, y) and applying the domain integration, we can obtain a weak-form integral equation: 𝑏 𝑎
𝑏 𝑎
∫ ∫ 𝜀𝑥 (𝑥, 𝑦)𝑃 (𝑥, 𝑦)𝑑 𝑥𝑑 𝑦 = ∫ ∫ 0 0
0 0
𝜕𝑈 (𝑥, 𝑦) 𝑃 (𝑥, 𝑦)𝑑 𝑥𝑑 𝑦 𝜕𝑥
Dividing the range (0, b) into 2n equal intervals. The interval length is: ℎ=𝑏∕2𝑛 And we denote: 𝑘𝜋 ℎ 𝑏
(4)
=𝜃 𝑟ℎ = 𝑦𝑟 𝑈 (𝑎, 𝑟ℎ) = 𝑈𝑟
(5)
where r is an integer. After a series of derivations, we have: ( ) [ ( ) ] 𝑏 𝑘𝜋𝑦 𝑘𝜋𝑏 + 𝛽𝑆2𝑟 + 𝛾𝑆2𝑟−1 𝑈 (𝑎, 𝑦)𝑑𝑦 = ℎ 𝛼𝑈 (𝑎, 𝑏) sin ∫ cos 𝑏 𝑏 0
Integrating the right side of Eq. (4) by parts we have: ∫0𝑏 ∫0𝑎 𝜀𝑥 (𝑥, 𝑦)𝑃 (𝑥, 𝑦)𝑑 𝑥𝑑 𝑦 = (𝑥,𝑦) 𝑑𝑥𝑑𝑦 ∫0𝑏 [𝑃 (𝑎, 𝑦)𝑈 (𝑎, 𝑦) − 𝑃 (0, 𝑦)𝑈 (0, 𝑦)]𝑑𝑦 − ∫0𝑏 ∫0𝑎 𝑈 (𝑥, 𝑦) 𝜕𝑃𝜕𝑥
It should be noted that, through integration by parts, the necessity of differentiating U is avoided, in the price of differentiating the test function P. In this way, we only need discrete values U when evaluating Eq. (4), so that the problem of numerical differentiation is transferred to numerical integration. In this work, ɛx (x, y) is expanded in Fourier series. Through expanding the region of interest to {(𝑥, 𝑦)|−𝑎 ≤ 𝑥 ≤ 𝑎, −𝑏 ≤ 𝑦 ≤ 𝑏} and imposing symmetry of ɛx (x, y) (i.e.𝜀𝑥 (𝑥, 𝑦) = 𝜀𝑥 (−𝑥, 𝑦) = 𝜀𝑥 (𝑥, −𝑦) = 𝜀𝑥 (−𝑥, −𝑦)), we can simplify the complete Fourier series as the cosine series: 𝜀𝑥 (𝑥, 𝑦) =
𝑁 ∑ 𝑀 ∑ 𝑛=0 𝑚=0
ℎ𝑛𝑚 cos
(
) ( 𝑚𝜋𝑦 ) 𝑛𝜋𝑥 cos 𝑎 𝑏
𝑆2𝑟 =
𝑟=0
( ) ( ) 𝑘𝜋𝑦2𝑟 𝑈 𝑎, 𝑦2𝑟 cos − 12 𝑈 (𝑎, 0) − 12 𝑈 (𝑎, 𝑏) cos 𝑘𝜋 𝑏
𝑛 ∑
𝑟=1
( ) 𝑈 𝑎, 𝑦2𝑟−1 cos
𝑘𝜋𝑦2𝑟−1 𝑏
(14)
when k is not equal to 0: 𝛼= 1∕𝜃+cos𝜃sin𝜃∕𝜃 2 − 2sin2 𝜃∕𝜃 3 𝛽= 2[(1 + cos2 𝜃)∕𝜃 2 − 2sin𝜃cos𝜃∕𝜃 3 ( ) 𝛾= 4 sin𝜃∕𝜃 3 − cos𝜃∕𝜃 2
(6)
when k is equal to 0: 𝛼= 0 𝛽= 2∕3 𝛾= 4∕3 Details about Filon’s integration can be found in [30]. 4. Verification using numerical experiments To investigate the performance of the weak-form method proposed in this paper, five cases of numerical experiments were conducted in this study: (1) a small homogeneous deformation, (2) two sinusoidal displacement fields, (3) a Gaussian displacement field deformations, (4) an analogous sinusoidal-Gaussian displacement field and (5) deformation in a perforated plate. In the case (1–4), computer-simulated images are used and the algorithm of generating computer-simulated speckle images is given in [31,32]. The simulated images adopt the following parameters: the size of simulated speckle images is 601 × 601 pixels, the number of speckle granules is 13,000 and the size of speckles is 2.5 pixels. Gaussian noises are added to both reference and deformed images to simulate the image noises encountered in practical experiments. The displacements, which are at regularly distributed 17,161 (131 × 131) points with a grid step of 4 pixels for the case (1) and 75,625 (275 × 275) points with a grid step of 2 pixels for the case (2–4), are computed using the IC–GN algorithm implemented in the software Ncorr taken from the website www.ncorr.com. In the case (5), the noisy displacement is created by adding normally distributed random numbers to the analytical displacement field. Then the proposed weak-form method is used to calculate the strain fields. In order to compare the accuracy of the proposed weak-form method with the LH method, we also use the LH method with C2 elements to calculate the strain fields of cases (3–5). All parameters of the LH method adopted in this paper can be found in [27]. To quantitatively investigate the performance of the proposed weakform method, the root-mean-square (RMS) error is calculated. The RMS error is defined as follows:
(8)
Substituting Eq. (6) into Eq. (8), and by using the orthogonality: ( ) ( ) 𝑗𝜋𝑦 𝑘𝜋𝑦 𝑏 cos 𝑑𝑦 = 𝛿𝑗𝑘 (9) ∫ cos 𝑏 𝑏 2 0
[ ( ) ( ) ] 𝑒𝑛𝑚 = ∫0𝑏 cos (𝑛𝜋) cos 𝑚𝜋𝑦 𝑈 (𝑎, 𝑦) − cos 𝑚𝜋𝑦 𝑈 (0, 𝑦) 𝑑𝑦+ 𝑏) 𝑏 ( ( ) sin 𝑛𝜋𝑥 cos 𝑚𝜋𝑦 𝑑 𝑥𝑑 𝑦 ∫0𝑏 ∫0𝑎 𝑈 (𝑥, 𝑦) 𝑛𝜋 𝑎 𝑎 𝑏
𝑛 ∑
𝑆2𝑟−1 =
𝑏
we have: 4𝜆 𝑛 ℎ𝑛𝑚 = 𝑒 𝑎𝑏 𝑚
(13)
where
where ℎ𝑛𝑚 are unknown coefficients. And (𝑁 + 1) × (𝑀 + 1) is called the number of expanding terms of the cosine series in this paper. According to the numerical results and our experience, we suggest that the range 𝐼 𝐽 of expanding terms is 15 × 15 ≤ (𝑁 + 1) × (𝑀 + 1) ≤ 𝐼3 × 𝐽3 , where I × J is the number of points calculated in the region of interest. To get a closedform expression of ℎ𝑛𝑚 , we utilize the orthogonality of cosine functions and choose test function P(x, y) as follows: ( ) ( ) 𝑗𝜋𝑥 𝑘𝜋𝑦 𝑃 (𝑥, 𝑦) = cos cos (7) 𝑎 𝑏 Substituting P(x, y) into Eq. (5), we have: ( ) ( ) cos 𝑘𝜋𝑦 𝑑 𝑥𝑑 𝑦 ∫0𝑏 ∫0𝑎 𝜀𝑥 (𝑥, 𝑦) cos 𝑗𝜋𝑥 𝑎 𝑏 [ ( ) ( ) ] 𝑘𝜋𝑦 𝑏 = ∫0 cos (𝑗𝜋) cos 𝑏 𝑈 (𝑎, 𝑦) − cos 𝑘𝜋𝑦 𝑈 (0, 𝑦) 𝑑𝑦+ 𝑏 ( ) ( ) cos 𝑘𝜋𝑦 𝑑 𝑥𝑑 𝑦 sin 𝑗𝜋𝑥 ∫0𝑏 ∫0𝑎 𝑈 (𝑥, 𝑦) 𝑗𝜋 𝑎 𝑎 𝑏
(12)
(10)
(11)
where, 𝜆 = 1∕4 when 𝑚 = 0 and 𝑛 = 0; 𝜆 = 1∕2 when one of m and n is 0; 𝜆 = 1 when both of m and n is not equal to 0. Strains distributions can be calculated by directly evaluating the integral of Eq. (11). However, traditional numerical integration algorithms, like the trapezoidal rule and the Gaussian integration, will give poor results, due to the sinusoidal terms in the integrand. Thus, Filon’s integration method is used which are specifically designed for integrands with an oscillatory factor of sin (kx)orcos (kx). 497
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 2. The x-directional displacement filed obtained by the IC–GN algorithm for the case of homogeneous deformation.
√ 𝑅𝑀𝑆 =
∑
(𝜀𝑚 − 𝜀)2 𝑛
(15)
where ɛm denotes the strain calculated by the PLS algorithm, or the proposed weak-form method at each calculation point; ɛ denotes the exact value of strains; n is the number of calculation points. 4.1. Homogeneous deformation As the simplest case of deformations with small strain-gradients, a homogeneous deformation of uniaxial tension is applied in x-direction: 𝑈 = 0.002𝑥, 𝜀𝑥 = 0.002, 0 ≤ 𝑥 ≤ 600, while the displacement V in y-direction is set to 0. Gaussian noise with a mean value of zero and a standard deviation of 0.5% of the entire 256 grayscales is added to both reference and deformed images. Fig. 2 shows the displacement field U obtained by the IC–GN algorithm. To obtain a reliable strain field, the proposed weak-form method is used to post-process the displacement calculated by the IC–GN algorithm. Fig. 3(a) gives the strain field obtained by the weak-form method with 9 × 9expanding terms (the RMS error is 98.36𝜇ɛ). For comparison, the strain field calculated using the PLS algorithm with the best strain calculation window of 21 × 21points (the RMS error is 70.93𝜇ɛ) is also shown in Fig. 3(b). Further, the effects of expanding terms of the proposed weakform method and the calculation windows size of the PLS algorithm are investigated. Fig. 4 shows the RMS errors of strain fields calculated by the proposed weak-form method with expanding terms ranging from 9 × 9 to 18 × 18, and by the PLS algorithm with various strain calculation windows varying from 3 × 3 points to 21 × 21 points. It can be found that the error dispersion of the weak-form method with 10 different expanding terms is smaller than that of the PLS algorithm with 10 different strain calculation windows. For the homogeneous deformation, the small number of expanding terms for the weak-form method will give a more accurate result, while large strain calculation windows are good choices for the PLS algorithm.
Fig. 3. The x-directional strain field obtained by (a) the weak-form method, and (b) the PLS algorithm, for the case of homogeneous deformation.
Example 1: The imposed displacement field of the first example is as follows: 𝑈 = sin
(
) ( ) 2𝜋 2𝜋𝑥 2𝜋𝑥 , 𝜀𝑥 = , 0 ≤ 𝑥 ≤ 600, cos 300 300 300
Gaussian noise with a mean value of zero and a standard deviation of 1% of the entire 256 grayscales is also added to both reference and deformed images. Fig. 5 shows the displacement field obtained by the IC–GN algorithm. Fig. 6 shows the strain field calculated by different algorithms. The side view of the strain field obtained by using the weak-form method with18 × 18 expanding terms (the RMS error is629.69𝜇ɛ) is shown in Fig. 6(a). Also, the side view of the strain field calculated using the PLS algorithm with the best strain window of 11 × 11 points (the RMS error is758.26𝜇ɛ) is shown in Fig. 6(b). For comparison, the side view of the exact strain field is shown in Fig. 6(c). Further, the effects of expanding terms of the proposed weak-form method and the strain calculation windows size of the PLS algorithm are investigated. Fig. 7 shows the RMS errors of strain fields calculated by the proposed weak-form method with expanding terms ranging from 18 × 18 to 27 × 27, and by the PLS algorithm with various strain calculation windows varying from 3 × 3 points to 21 × 21 points.
4.2. Sinusoidal displacement fields To study the performance of the proposed weak-form method for large strain-gradient or highly-oscillatory deformations, the next two examples are considering in-plane sinusoidal displacement fields in the x-direction. 498
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 4. RMS errors of the estimated strain fields with various expanding terms of the weak-form method and various strain calculation windows of the PLS algorithm, for the case of homogeneous deformation.
Fig. 5. The x-directional displacement field obtained by the IC–GN algorithm for example 1 of the sinusoidal displacement fields.
Example 2: The second example with larger strain gradients and larger noise than the first example is given as follows: ( ) ( ) 2𝜋 2𝜋𝑥 2𝜋𝑥 𝑈 = sin cos , 𝜀𝑥 = , 0 ≤ 𝑥 ≤ 600, 200 200 200 Gaussian noise with a mean value of zero and a standard deviation of 2% of the entire 256 grayscales is also added to both reference and deformed images. Fig. 8 shows the computed x-direction displacement field using the IC–GN algorithm. By post-processing the displacement calculated by the IC–GN algorithm using the weak-form method with 18 × 18 expanding terms, the side view of the strain field is shown in Fig. 9(a). Besides, the side view of the strain field calculated using the PLS algorithm with the best strain window of 11 × 11 points is shown in Fig. 9(b). For comparison, the side view of the exact strain field is shown in Fig. 9(c). In a similar manner, the errors of strain fields corresponding to various expanding terms of the weak-form method and different strain calculation windows of the PLS algorithm are shown in Fig. 10. From these two examples, it is seen that the proposed method has advantages for deformations with large strain-gradients. The first ad-
Fig. 6. Side views of strain fields obtained by (a) the weak-form method, (b) the PLS algorithm, and (c) the exact strain, for example 1 of the sinusoidal displacement fields. 499
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 7. RMS errors of the estimated strain fields with various expanding terms of the weak-form method and various strain calculation windows of the PLS algorithm, for example 1 of the sinusoidal displacement fields.
Fig. 8. The x-directional displacement field obtained by the IC–GN algorithm for example 2 of the sinusoidal displacement fields.
vantage of the weak-form method is that the error dispersion of strains with 10 different expanding terms is smaller than the PLS algorithm with 10 different strain calculation windows. The other advantage is that strain fields calculated by the weak-form method with various expanding terms are more accurate and less insensitive to noise. 4.3. A Gaussian displacement field A Gaussian displacement field proposed by Xu et al. [33] is considered here: ( ) √ 𝑈 = 𝑎 × 𝑐 𝜋∕2𝑒𝑟𝑓 √𝑥 , −300 ≤ 𝑥 ≤ 300 2𝑐
𝜀𝑥 = 𝑎
2 −𝑥 × 𝑒 2𝑐 2
where parameter a is the intensity of the Gauss function and c is the standard error of the Gauss function. In this paper, a is 0.0025 and c is 16. Gaussian noise with a mean value of zero and a standard deviation of 1% of the entire 256 grayscales is added to both reference and deformed images.
Fig. 9. Side views of strain fields obtained by (a) the weak-form method, (b) the PLS algorithm, (c) the exact strain, for example 2 of the sinusoidal displacement fields. 500
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 10. RMS errors of the estimated strain fields with various expanding terms of the weak-form method and various strain calculation windows of the PLS algorithm, for example 2 of the sinusoidal displacement fields.
Fig. 11. The side view of the x-directional displacement field obtained by the IC–GN algorithm for the Gaussian displacement field.
Fig. 11 shows the computed x-direction displacement field using the IC–GN algorithm. By post-processing the obtained displacement using the weak-form method with 36 × 36 expanding terms, the side view of the strain field is shown in Fig. 12(a). Besides, the side views of the strain field calculated using the PLS algorithm with the best strain window of 7 × 7 points and using the LH method with the best strain window of 29 × 29 points are shown in Fig. 12(b) and Fig. 12(c), respectively. The errors of strain fields corresponding to various expanding terms of the weak-form method and different strain calculation windows of the PLS algorithm as well as the LH method are shown in Fig. 13. 4.4. An analogous sinusoidal-Gaussian displacement field An analogous sinusoidal-Gaussian displacement field proposed by Yang et al. [34] is considered here: 𝑈 (𝑥, 𝑦) = sin (2𝜋𝐸 (𝑥)) sin (2𝜋𝐸 (𝑦))∕10, ) 2 ( 𝐸 (𝑥) = 𝑒−(𝑥−300) ∕ 2𝜎 2 ) 2 ( − ( 𝑦 −300) 2 𝐸 (𝑦 ) = 𝑒 ∕ 2𝜎
0 ≤ 𝑥 ≤ 600 Fig. 12. Side views of strain fields obtained by (a) the weak-form method, (b) the PLS algorithm, (c) the LH method, and (d) the exact strain, for the Gaussian displacement field. 501
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 13. RMS errors of the estimated strain fields with various expanding terms of the weak-form method, various strain calculation windows of the PLS algorithm and the LH method, for the Gaussian displacement field.
where 𝜎 denotes the Gaussian Root-Mean-Square (RMS) width, which is 80 in this paper. Gaussian noise with a mean value of zero and a standard deviation of 1% is added to both reference and deformed images. Fig. 14 shows the side views of the computed x-direction displacement field using the IC–GN algorithm. By post-processing the displacement using the weak-form method with 27 × 27 expanding terms, the side view of the strain field is shown in Fig. 15(a). Besides, the side views of strain field calculated using the PLS algorithm with the best strain window of 7 × 7 points and using the LH method with the best strain window of 29 × 29 points are shown in Fig. 15(b) and Fig. 15(c), respectively. The errors of strain fields corresponding to a various number of expanding terms of the weak-form method and different strain calculation windows of the PLS algorithm as well as the LH method are shown in Fig. 16.
4.5. Deformation in a perforated plate To validate the feasibility of the proposed method for perforated plate problems, the last case investigated is a plate with a small hole under a uniform tension in x-direction, as shown in Fig. 17. The material properties of the plate are 𝐸 = 0.8𝐺𝑃 𝑎 and 𝑣 = 0.3. The radius of the hole is 50mm and the tension q is 23 × 107 𝑁∕𝑚2 . The size of the region of interest is 598 × 598mm2 with a grid step of 2mm. The analytical solutions of the displacement field and the strain field can be found in [35], and the noisy displacement is created by adding normally distributed random numbers to the analytical displacement field. In this example, the weak-form method described in Section 3.1 cannot be directly used, because of the existence of the hole in the center of the plate. A simple pre-processing of the displacement data is implemented before calculating the strain field: dividing the original image into 4 subdomains, and extrapolating displacements in each subdomain to a rectangular box. Fig. 18(a) illustrates an irregular region of interest, in which displacements at red points are already calculated by NR, IC–GN or other algorithms. In Fig. 18(b), the virtual displacements at blue points are to be extrapolated by using displacements at red points. The dataextrapolating process is described as follows: (1) a square calculation window containing 5 × 5 − 1 red points and 1 blue point (point P) is chosen to calculate the virtual displacement at the point P, as shown in Fig. 18(b). (2) The displacement field in the calculation window can
Fig. 14. (a)Theoretical x-directional displacement in 3D view, (b) theoretical x-directional displacement in side view, and (c) x-directional displacements obtained by the IC–GN algorithm in side view.
502
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 16. RMS errors of the estimated strain fields with various expanding terms of the weak-form method, various strain calculation windows of the PLS algorithm and the LH method, for the analogous sinusoidal-Gaussian displacement field.
Fig. 17. A plate with a small hole submitted to a uniform tension in x-direction.
be approximated by using the polynomial 𝑈 = 𝑎0 + 𝑎1 𝑥 + 𝑎2 𝑦 + 𝑎3 𝑥2 + 𝑎4 𝑥𝑦 + 𝑎5 𝑦2 , in which the unknown polynomial coefficientsa0 , a1 , a2 , a3 , a4 , a5 are solved by the least-squares method using displacements at blue points. (3) After the virtual displacement at the point P is obtained, it will be treated as a known value to calculate the virtual displacement at the next blue point P1 . The same procedure is used to calculate virtual displacements at each blue point from left to right and from top to bottom in turn. Finally, when all virtual displacements at blue points are calculated, the proposed weak-form method can be used to calculate the strain field in the rectangular region shown in Fig. 18(b). For the perforated plate problem, it should be noted that the order of polynomial using in extrapolating process will affect the accuracy of the strain to some extent. By dividing the region of interest into 4 subdomains as shown in Fig. 17, the strain field can be calculated by using the process described above in each subdomain. Fig. 19 shows ɛx calculated using the proposed weak-form method with 30 × 30expanding terms in each subdomain. Fig. 20 shows the ɛx along the y-axis direction calculated by the proposed weak-form method, the PLS algorithm, and the analytical solution, respectively. Fig. 21 displays the RMS error of strain field corresponding to various expanding terms of the weak-form method, different strain calculation windows of the PLS algorithm and the LH method. It is observed that the weak-form method can improve the accuracy of
Fig. 15. Side views of strain fields obtained by (a) the weak-form method, (b) the PLS algorithm, and (c) the LH method for the analogous sinusoidal-Gaussian displacement field.
503
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Fig. 19. The x-directional strain field of the perorated plate submitted to a uniform tension obtained by the weak-form method.
Fig. 20. The x-directional strain field along the y-axis of the perorated plate obtained by the weak-form method, the PLS algorithm with the strain window of 5 × 5 as well as 11 × 11 points, and the exact value.
Fig. 18. (a) an irregular region of interest, (b) a rectangle region by extrapolating displacements for the irregular region of interest. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)
strain fields especially in obvious stress concentrated areas compared to the PLS algorithm. In order to elucidate the computational efficiency of the proposed method, the computational time of above all numerical experiments on the condition of the best-expanding terms and the best strain window is shown in Table 1. It is seen that the proposed method has similar computational efficiency as compared to the PLS algorithm in most of these numerical simulations.
Fig. 21. RMS errors of the estimated strain fields with various expanding terms of the weak-form method, various strain calculation windows of the PLS algorithm and the LH method, for the perorated plate. 504
L. Dong, Q. Yang and J. Wang
Optics and Lasers in Engineering 121 (2019) 495–505
Table 1 The computational time for all numerical experiments on the condition of the best expanding terms of the weak-form method and the best strain windows of the PLS algorithm. EF is the expanding terms of the weak-form method. SCW is the strain calculation windows of the PLS algorithm. Method time (ET or SCW) Numerical experiments
The weak-form method
The PLS algorithm
Homogeneous deformation Example 1 of the sinusoidal displacement field Example 2 of the sinusoidal displacement field The Gaussian displacement field The analogous sinusoidal-Gaussian displacement field The deformation in a perforated plate
0.33 s (ET:9 × 9) 4.02 s (ET:18 × 18) 4.23 s (ET:18 × 18) 10.71 s (ET:29 × 29) 9.10 s (ET:27 × 27) 11.04 s (ET:30 × 30)
14.31 s (SCW:21 × 21) 20.29 s (SCW:11 × 11) 20.71 s (SCW:11 × 11) 11.34 s (SCW:7 × 7) 11.23 s (SCW:7 × 7) 8.00 s (SCW:5 × 5)
5. Conclusion
[7] Bruck HA, McNeill SR, Sutton MA, Peters WH. Digital image correlation using Newton–Raphson method of partial differential correction. Exp Mech 1989;29:261–7. [8] Lu H, Cary PD. Deformation measurements by digital image correlation: implementation of a second-order displacement gradient. Exp Mech 2000;40:393–400. [9] Yoneyama S, Morimoto Y. Accurate displacement measurement by correlation of colored random patterns. JSME Int J Series A 2003;46:178–84. [10] Baker S, Matthews I. Lucas–Kanade 20 years on: a unifying framework. Int J Comput Vis 2004;56:221–55. [11] Wang CCB, Deng JM, Ateshian GA, Hung CT. An automated approach for direct measurement of two-dimensional strain distributions within articular cartilage under unconfined compression. J Biomech Eng 2002;124:557–67. [12] Sutton MA, Turner JL, Bruck HA, Chae TA. Full-field representation of discretely sampled surface deformation for displacement and strain analysis. Exp Mech 1991;31:168–77. [13] Tong W. Detection of plastic deformation patterns in a binary aluminum alloy. Exp Mech 1997;37:452–9. [14] Yoneyama S. Smoothing measured displacements and computing strains utilising finite element method. Strain 2011;47:258–66. [15] Han Y, Kim DW, Kwon HJ. Application of digital image cross-correlation and smoothing function to the diagnosis of breast cancer. J Mech Behav Biomed Mater 2012;14:7–18. [16] Grama SN, Subramanian SJ. Computation of full-field strains using principal component analysis. Exp Mech 2014;54:913–33. [17] Zhao J, Song Y, Wu X. Fast Hermite element method for smoothing and differentiating noisy displacement field in digital image correlation. Opt Lasers Eng 2015;68:25–34. [18] Wang D, DiazDelaO FA, Wang W, Lin X, Patterson EA, Mottershead JE. Uncertainty quantification in DIC with Kriging regression. Opt Lasers Eng 2016;78:182–95. [19] Wattrisse B, Chrysochoos A, Muracciole J-M, Némoz-Gaillard M. Analysis of strain localization during tensile tests by digital image correlation. Exp Mech 2001;41:29–39. [20] Pan B, Xie H, Guo Z, Hua T. Full-field strain measurement using a two-dimensional Savitzky–Golay digital differentiator in digital image correlation. Opt Eng 2007;46:033601. [21] Pan B, Asundi A, Xie H, Gao J. Digital image correlation using iterative least squares and pointwise least squares for displacement field and strain field measurements. Opt Lasers Eng 2009;47:865–74. [22] Lehoucq RB, Reu PL, Turner DZ. A novel class of strain measures for digital image correlation. Strain 2015;51:265–75. [23] Dai X, Yang F, Chen Z, Shao X, He X. Strain field estimation based on digital image correlation and radial basis function. Opt Lasers Eng 2015;65:64–72. [24] Pan B, Yuan J, Xia Y. Strain field denoising for digital image correlation using a regularized cost-function. Opt Lasers Eng 2015;65:9–17. [25] Schreier HW, Sutton MA. Systematic errors in digital image correlation due to undermatched subset shape functions. Exp Mech 2002;42:303–10. [26] Hwang SF, Wu WJ. Deformation measurement around a high strain-gradient region using a digital image correlation method. J Mech Sci Technol 2012;26:3169–75. [27] Li X, Fang G, Zhao JQ, Zhang ZM, Wu XX. Local Hermite (LH) Method: an accurate and robust smooth technique for high-gradient strain reconstruction in digital image correlation. Opt Lasers Eng 2019;112:26–38. [28] Liu CS, Dong L. Closed-form higher-order numerical differentiators for differentiating noisy signals. Appl Math Comput 2019 in press. [29] Dong L, Alotaibi A, Mohiuddine SA, Atluri SN. Computational methods in engineering: a variety of primal & mixed methods, with global & local interpolations, for well-posed or ill-posed BCs. Comput Model Eng Sci 2014;99:1–85. [30] Filon LNG III. On a quadrature formula for trigonometric integrals. In: Proceedings of the Royal Society of Edinburgh, 49; 1930. p. 38–47. [31] Zhou P, Goodson KE. Subpixel displacement and deformation gradient measurement using digital image/speckle correlation (DISC). Opt Eng 2001;40:1613–20. [32] Pan B, Xie HM, Xu BQ, Dai FL. Performance of sub-pixel registration algorithms in digital image correlation. Measur Sci Technol 2006;17:1615–21. [33] Xu X, Su Y, Cai Y, Cheng T, Zhang Q. Effects of various shape functions and subset size in local deformation measurements using DIC. Exp Mech 2015;55:1575–90. [34] Yang X, Chen X, Xi J. Comparative analysis of warp function for digital image correlation-based accurate single-shot 3D shape measurement. Sensors 2018;18:1208. [35] Timoshenko SP, Goodier JN. Theory of elasticity; 1970.
The assessment of strain fields is important for deformation measurement. However, the strain fields obtained directly by the NR algorithm, the IC–GN algorithm, or genetic algorithms, have poor accuracy. The PLS algorithm is widely used for retrieving strain distribution. However, the PLS algorithm has some limitations when dealing with large gradients strain fields, especially with noisy measurements. In this paper, an alternative strain estimation approach based on the weak form is proposed, and verified by numerical simulations. Numerical simulations reveal that when strain gradients are small, strain fields obtained by the weak-form method are comparable to the PLS algorithm. However, for large strain-gradient or highly-oscillating deformations, the results obtained by the weak-form method are more accurate than the PLS algorithm. Furthermore, the error dispersion of the weak-form method is smaller than the PLS algorithm no matter strain gradients are small or large. For complex large strain-gradient or highlyoscillating deformation problems, although the LH method can obtain more accurate strain fields than the PLS algorithm, it is much more computationally expensive and complex than the weak-form method and the PLS algorithm. Through numerical experiments, one can also see that the number of expanding terms affects the accuracy of the weak-form method, although in a certain range the error dispersion is quite small. In the 𝐼 𝐽 suggested range [ 15 × 15 , 𝐼3 × 𝐽3 ] for the number of expanding terms, a small number is preferred for small strain-gradients and for large noises. It should also be noted that, while the weak-form method described in Section 3.1 is directly used to calculate the strain distribution for a rectangular region, irregular regions of interest can also be considered by extrapolating the original displacement data to a rectangular box. This is demonstrated in Section 4.5 with the numerical simulation of a perforated plate. Acknowledgment The authors thankfully acknowledge the support from the National Key Research and Development Program of China (No. 2017YFA0207800). References [1] Nakadate S, Yatagai T, Saito H. Electronic speckle pattern interferometry using digital image processing techniques. Appl Opt 1980;19:1879–83. [2] Bing P, Kemao Q, Huimin X, Anand A. Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review. Measur Sci Technol 2009;20:062001. [3] Pan B. Recent progress in digital image correlation. Exp Mech 2011;51:1223–35. [4] Peters WH, Ranson WF. Digital imaging techniques in experimental stress analysis. SPIE; 1982. p. 5. [5] Peters WH, Ranson WF, Sutton MA, Chu TC, Anderson J. Application of digital correlation methods to rigid body mechanics. SPIE; 1983. p. 5. [6] Pan B. Reliability-guided digital image correlation for image deformation measurement. Appl Opt 2009;48:1535–42.
505