A hybrid DIC–EFG method for strain field characterization and stress intensity factor evaluation of a fatigue crack

Measurement 154 (2020) 107498

Contents lists available at ScienceDirect

Measurement journal homepage: www.elsevier.com/locate/measurement

A hybrid DIC–EFG method for strain field characterization and stress intensity factor evaluation of a fatigue crack Zhu Zhihui a,b,⇑, Luo Sihui a, Feng Qianshuo a, Chen Yuanchang c, Wang Fan a, Jiang Lizhong a a

School of Civil Engineering, Central South University, Changsha 410075, China National Engineering Laboratory for High Speed Railway Construction, Central South University, Changsha 410075, China c Structural Dynamics and Acoustic Systems Laboratory, University of Massachusetts Lowell, Lowell, MA 01854, USA b

a r t i c l e

i n f o

Article history: Received 28 July 2019 Received in revised form 6 December 2019 Accepted 4 January 2020 Available online 13 January 2020 Keywords: Digital image correlation Element-free Galerkin Meshless Strain field Stress intensity factor J-integral Fracture performance Fatigue crack

a b s t r a c t A Hybrid DIC–EFG Method combining Digital Image Correlation (DIC) and Element-free Galerkin (EFG) method is proposed to characterize the strain field and study the fracture performance of a crack by evaluating the stress intensity factor (SIF). The full-field displacement data is measured by DIC. The meshless calculation model is developed by using the speckle image coordinate. Then the displacement data is imported into the model and strain is calculated by the EFG method. The SIF is obtained based on the J-integral method. A static tensile test with a quasi-static cyclic load is conducted on a Compact Tensile specimen (CT specimen) to experimentally validate the accuracy and stability of proposed method. The proposed method can effectively predict the strain at the crack tip, which is a great challenge for the DIC measurement. Compared with the reference ASTM standard, the proposed method can obtain the SIF accurately without knowing the geometry, load and boundary condition of the structure. Ó 2020 Elsevier Ltd. All rights reserved.

1. Introduction In civil and bridge applications, some structures have irregular gaps, such as the U-shaped notch of the diaphragm in the orthotropic plate. Under the continuous traffic loads and excitation, the tip of the notch tends to have fatigue issues, and the cracks due to the fatigue would sometimes generate [1,2]. These fatigue cracks can decrease the service life of the components and may lead to structural fracture or damage [3]. Thus, it is important to monitor the structural health and accurately evaluate the fracture performance of the fatigue cracks for the safety and durability of structures. The Stress Intensity Factor (SIF) is used to characterize the stress state near the crack tip under the external load or residual stresses. It is an important indicator of fracture performance [4]. The SIF can be obtained by analytical method, numerical simulation or experimental method. The analytical method can be only applied to the structures with simple boundaries and regular geometry. For the cracks with geometric discontinuity and nonlinear material, the numerical simulation method can be used, but it has some limitations that the shape and the geometry of the crack ⇑ Corresponding author at: School of Civil Engineering, Central South University, Changsha 410075, China. E-mail address: [email protected] (Z. Zhu). https://doi.org/10.1016/j.measurement.2020.107498 0263-2241/Ó 2020 Elsevier Ltd. All rights reserved.

are difficult to model for the numerical simulation [5]. The experimental methods for SIF have received much attention, including the strain gauge method [6] and the light measurement method [7,8]. The strain gauge method has strict requirements on the placement and angle of the patch, and the patch is complicated with limited measurement points; the compliance method cannot be used for dynamic measurement; and the optical method has been applied to the calculation of fracture mechanics SIF because of its characteristics of full-field measurement information and dynamic measurement [9,10]. Roux et al. [7] estimated the SIF of the two-dimensional and three-dimensional crack based on digital image interference technology. Hadadian et al. [11] calculated the SIF of mixed I-II crack based on electronic speckle interferometry. The traditional full-field displacement and strain photometric methods include photoelastic method, moire interferometry method and electronic speckle pattern interferometry [12], etc. These methods require pretreatment of the test piece surface and a precise optical path. The electronic speckle pattern interferometry requires strict darkroom and vibration isolation conditions [12,13]. The Digital Image Correlation (DIC) measures the surface displacement of the test piece by comparing the grayscale of the original image with the target image [14,15], which has low requirements for the test environment. The measurement set up is simple without complicated optical imaging systems [16]. In

2

Z. Zhu et al. / Measurement 154 (2020) 107498

recent years, DIC has been widely used in structural dynamics, structural health monitoring, and fracture mechanics. Shah et al. [17] experimentally validated that DIC was a very useful tool for determining the fracture parameters such as crack opening displacement, crack length and crack tip position. Fawad et al. [18] proved that the SIF obtained from the DIC had a good agreement with the finite element numerical calculation as well as the theoretical solution. Chen et al. [19,20] proposed an expansion procedure to construct the full-field displacement and strain information by using either a model approach or a non-model approach. The methods for calculating the SIF based on DIC include the least square fitting method [21], J-integral method [22] and the DIC-FEM coupling method [5]. Some analyses used least-squares methods to fit the Williams’ series to the DIC results and obtained the mode I SIF [23]. However, Roux and Bhardwaj et al. [24,25] stated that the fitting method was sensitive to the position of the crack tip, analyzing critical locations can become challenging. The J-integral method holds for any integration contour surrounding the crack tip, either linear elastic or elastoplastic. Hence, its formulation does not require the accurate location of the crack tip, which is a significant source of uncertainty for other approaches [26]. But the DIC technique is difficult to obtain the measurement results around the crack and the edge, and some data will be lost when the SIF is calculated based on the J-integral. Barhli [14] presented that the information of a crack tip could be effectively obtained by coupling DIC-FEM to SIF. But due to the crack tip stress singularity, it often requires an extensive number of degrees of freedom, meshing is troublesome. The Element Free Galerkin (EFG) method has the characteristics of high accuracy and fast convergence compared with FEM [27]. It can also be used to create a meshless model by using the DIC measured results. The preprocessing procedure is simple. The meshless model does not need to re-mesh the crack tip, which does not have the grid distortion issue [28]. Therefore, a full-plane model including crack information can be established by coupling the DIC and EFG methods, and the accurate full-plane strain field can be calculated. The hybrid DIC–EFG method can obtain the accurate full-field strain including the crack tip, which is a great challenge for the DIC. And then SIF is calculated by the J-integral method. Compared with the east square fitting method, the hybrid DIC-EFG is more stable at the crack tip. This paper proposes a hybrid DIC–EFG method based on DIC and EFG for extracting the SIF of a surface crack. The DIC is used to obtain the full-field displacement and node information and then develop a meshless model. The EFG method can directly obtain full-field stress and strain. Then the J-integral method is used to calculate the SIF of a fatigue crack. Based on the VIC-3D instrument, a quasi-static cyclic static tensile test is performed on a compact tensile specimen (CT) with Q235 material to validate the method experimentally. The sensitive analysis were described to studed the parameter value and stability of proposed method.

2. SIF calculation method based on the EFG method

Fig. 1. Theoretical basis of the meshless method.

domain needs to be approximated by using the function value of the field node in the locally supported domain. Let uh ðxÞ be an approximate solution, which is called the trial-function:

uðxÞ  uh ðxÞ ¼

m X

pi ðxÞai ðxÞ  pT ðxÞaðxÞ

where pi ðxÞ is the basis function, m is the number of the term of pi ðxÞ, aðxÞ is a specific coefficient matrix, which is related to the spatial coordinate x. The basic function pðxÞ can be selected according to the two-dimensional Pascal triangle. The larger m is, the more accurate the trial function is Fig. 2. The function is solved based on the moving least squares method (MLS). In the case of the smallest error, the extreme value of the trial function is required, and the approximate displacement solution of the EFG method is obtained as:

uh ðxÞ ¼

nN X

NkI ðxÞuI  N k u

ð2Þ

I¼1

where N k is a shape function of the EFG method, and the specific expression is:

N k ¼ pT ðxÞA1 ðxÞBðxÞ A ¼ P T W ðxÞP; B ¼ P T WðxÞ

ð3Þ

It can be seen that the shape function is mainly affected by the weight function W ðxÞ and the basis function pðxÞ. Appropriate weight function and basis function need to be selected to solve the problem. The spline function is often selected as a weight function to take values. The cubic spline function is:

W i ðxÞ ¼

8 > <

2 3

 4r 2i þ 4r 3i r i  0:5

4  4r i þ 4r 2i  43 r 3i 0:5 < r i  1 > :3 0r i > 1

ð4Þ

where r i ¼ s0diðxÞ, s0 ðxÞ ¼ kx  xI k is the distance between the field node x and the sampling point xI . The definition of each distance is as shown in Fig. 3 and di is the search radius of the EFG method, i.e., the weight function Support domain size. Considering the particularity of the crack tip, when calculating the strain, the node search for the shape function should be

2.1. EFG method The meshless method is a method to establish algebraic system equations based on nodes without a grid. The calculation model is represented by a set of nodes scattered over it. Because these nodes are the carriers of meshless field variables, they are called field nodes. The node density can be adjusted according to the requirement of calculation accuracy Fig. 1. As there is no grid in the meshless method, the variable u (displacement, stress, etc.) at arbitrary point X= (x, y, z) in the problem

ð1Þ

i¼1

Fig. 2. 2D Pascal Triangle.

3

Z. Zhu et al. / Measurement 154 (2020) 107498

calculation, the path of the J-integral is selected as a rectangle centered on the crack. J-integral [29] are defined as:

Z

J¼

Wdy  T C

@u ds @x

ð9Þ

where C is an arbitrary integral loop around the crack tip from one point counterclockwise to the upper surface of the crack, W is the strain energy density function, u is the displacement vector, and T is the stress tensor on the integral loop, and the integral The direction of the line is related to:

T¼

Fig. 3. Diffraction method.

searched by the diffraction method. The distance between the crack tip and the node around the crack to the integration point is:

sðxÞ ¼

 k s1 þ s2 ðxÞ s0 ðxÞ s0 ðxÞ

ð5Þ

among them, s1 ðxÞ ¼ kxc  xI k, s2 ðxÞ ¼ kx  xc k. The definition of each distance is as shown in Fig. 3 and k is a parameter, which is generally 1 or 2. The diffraction method eliminates discontinuities in the domain and is more suitable for non-convex boundaries. The displacement of the node is interpolated based on the MLS shape function of the EFG method (formula (3)) to, and the displacement result can be obtained according to the formula (2), and according to the mechanic’s principle The stresses and strains of the corresponding nodes can be obtained:

2

@ @x

e¼6 40

@ @y

0

3

@ 7 @y 5

 Nku

ð6Þ

@ @x

2 1 E 6 l r¼ 4 1  l2 0

l 1 0

0

3

0 7 5e

ð7Þ

1l 2

where e and r are the strain and stress vectors, respectively, N k is the shape function matrix of the field node corresponding to the EFG method, and E and l are the elastic moduli and Poisson’s ratio of the material. 2.2. SIF calculation method For standard CT specimens under static tensile loads, the calculation of the stress intensity factor yields an accurate analytical solution. According to the international ASTM standard, the stress intensity factor K I of the type I crack is:

P K I ¼ pffiffiffiffiffiffi FðaÞ B W F ðaÞ ¼ ð2 þ aÞ





rxx cos ðn; xÞ þ sxy cosðn; yÞ i   þ sxy cos ðn; xÞ þ ryy cosðn; yÞ j

where n is a straight line pointing to the normal vector inside the integral line, the path that can be selected is parallel to the coordinate axis, and the line integral is converted into the integral of the coordinates. For a two-dimensional uniform plane with I-type crack, the method for obtaining the stress intensity factor based on J-integration is:

Plane stress K I ¼

pffiffiffiffiffiffiffiffiffiffiffiffi JðCÞE

Plane stress K I ¼

ð11Þ

qffiffiffiffiffiffiffiffiffi JðCÞE 1l2

3. Hybrid DIC–EFG method DIC is an advanced optical measurement method based on digital image processing and numerical calculation. However, because DIC can’t determine the displacement vector at the discontinuity region of the displacement field, the measurement results are hard to capture the data at the crack tip. [14] The hybrid DIC–EFG method can obtain full-field information, including the crack tip. And the calculation of SIF is not sensitive to the position of the crack tip. 3.1. Digital image correlation (DIC) method In DIC, two digital images before and after the deformation are addressed. By searching the deformed image for the subset with the highest correlation of the original image subset, the surface displacement and strain of the test piece are calculated by comparing the grayscale.[14] The digital image before the deformation is generally referred to as ‘‘reference image”, and the deformed digital image is referred to as ‘‘post-deformed image”. Its basic principle is shown in Fig. 4. A rectangular reference image sub-region of size (2 m + 1) pixel(2 m + 1) pixel centered on the point (x, y) to be

ð8Þ

0:866 þ 4:64ðaÞ  13:32a2 þ 14:72a3  5:6a4 ð1  aÞ3=2

where a is the ratio a=W of the crack length a to the sample width W, and the function FðaÞ is the shape function of the boundary tension plane plate, where B is the thickness of the sample and P is the load size. The J-integral method is used to calculate the stress intensity factor of the model. The J-integral method integrates the encirclement line [26]. It is only necessary to know the stress and the displacement solution on the far-field integral path, and it can be integrated at any given line around the crack tip. Since the Jintegral is independent of the path, for the convenience of

ð10Þ

Fig. 4. Principle of digital image correlation method.

4

Z. Zhu et al. / Measurement 154 (2020) 107498

sought is taken in the reference image, and a deformed image is taken in the deformed image. In the sub-region, the deformed image sub-region and the reference image sub-region have a certain correlation, which can be expressed by correlation coefficients. According to a search method, a deformed image sub-region centered on (x’, y’) is found. The correlation reaches a maximum value, and the reference image sub-region corresponds to the deformed image sub-image after deformation. The displacement information of the point can be obtained according to the coordinate difference between the two points [30]. The DIC instrumentation VIC-3D is used here as an example. The correlation function C used is the sum of squared errors (SSD) of the pixel values in the reference image sub-region and the target image sub-region. The smaller the sum, the higher the similarity.

C ðx; y; u; v Þ ¼

n=2 X

ðIðx þ i; y þ jÞ  I ðx þ u þ i; y þ v þ jÞÞ

2

ð12Þ

i;j¼n=2

where x, y are the pixel coordinates in the reference image, u, v is the actual coordinate value, and Iðx þ i; y þ iÞ represents the gray value at the reference image ðx þ i; y þ jÞ, I ðx þ u þ i; y þ v þ jÞ represents the gray value at the target image ðx þ u þ i; y þ v þ jÞ. The size of the search sub-area recommended by the software depends on the image quality, including the grayscale distribution (the size and distribution of the scattered spots), the degree of clarity, the brightness, etc. The larger the subset, the more image information it contains, the more accurate the results will be, and the longer the calculation time. 3.2. Meshless computing model Fig. 5 shows the formation process of the meshless computational model S. The DIC analysis is used to retrieve the node distribution and displacement vector of the full-field range Uc. VIC-3D software is used for image recognition, and the least square coefficient (12) is used for correlation matching. The VIC-3D analysis software calculates the subset value of the correlation coefficient C according to the speckle quality. The analysis of step A is performed under the condition of selecting a specific subset size and

calculating steps. Searching with a large subset size is about 51  51 pixels to 81  81 pixels. As shown in Fig. 5(a), the selected steps are kept constant, finally, the displacement of the full-field U0 is obtained as shown in Fig. 5(b). Due to the presence of cracks, the discontinuity of the scattered spots disturbs the DIC analysis, thus losing the crack and the calculation results around the boundary of the sample. As shown in Fig. 5(c), the coordinate position (x, y) and the displacement (ux, uy) of the field node in the range S0 are derived based on the distribution of the field nodes at a certain pitch delta, which is the node spacing of the meshless model Fig. 6. The calculation result of the crack tip is performed using step B with a small subset, which is about 19  19 pixels, as shown in Fig. 5 (d). Small subsets can search for node distributions at boundary conditions and crack locations and calculate strain values, but provide less accurate displacement and strain fields. The cracked portion Sc displacement measurement result Uc is injected into the full-field S0 for filling the blank of the calculation result around the crack to form a displacement field of the full-field S, as shown in Fig. 5 (g). According to the encryption interval (0.2–0.5) *delta distribution field node, the coordinate position (x, y) and displacement (ux, uy) of the corresponding field node in the crack tip range Sc are derived. Therefore, step B is to obtain the crack tip displacement measurement result and field node arrangement information. Among them, ABCD is the supplementary part of the data points missing in the crack point of S0, supplemented by the Sc model, the node spacing is about delta/4, and A´ B´ C´ D´ is the transition area, which is composed of S0 area and Sc. The node spacing is approximately delta/2 to avoid abrupt stress changes due to displacement errors. The size of the transition area is related to the calculated search radius di of the EFG method, and the distance from the ABCD boundary should be not less than di. di is the search radius, which is determined by the calculation method of the EFG method and is generally 1 to 3 times the node spacing delta. di is independent of the transition area A’B’C’D’. The transition area A’B’C’D’ generally takes 1 to 2 times of the search radius di. Special attention should be paid to the coordinate system reference information in the sample before the measurement. The general situation is the original position and the X-axis direction. As far as possible, the crack is located on the X-axis. At the same time, the node information should be shared after each calculation.

Fig. 5. DIC Effective Displacement Field and Crack Meshless Computation Model.

5

Z. Zhu et al. / Measurement 154 (2020) 107498

Fig. 6. The process of the hybrid DIC–EFG method.

3.3. Calculation process of the hybrid DIC–EFG method The process of hybrid DIC–EFG method can be organized as follows: DIC retrieves the full-field information by using different subsets, developing the meshless calculation model and obtaining the full-field displacement field. EFG method is used to calculate the stress, strain field and extract the J-integral to calculate the SIF. The field nodes required for the meshless model can be derived from the DIC analysis. The nodes are distributed according to the distance delta  delta in the full-field range, and the node of the crack tip portion is refined appropriately. The DIC will automatically identify the crack tip position and derive the crack edge position. Finally, the meshless model S is formed by combining the fullfield area S0 field node information and the crack tip area Sc field node information. Based on the EFG method, combined with the displacement field Su , the strain field Se of S and the stress field Sr are calculated. The fracture field obtained from the measurements and theoretical analysis is based on the fracture behavior of the material using the EFG method. For the discrete point distribution node to calculate the J-integral, it is first necessary to obtain the displacement, stress, and strain of the discrete integral points on the integral line. Therefore, the MLS shape function based on the EFG method (formula (3)) is used to interpolate the discrete points on the integral line. For certain crack tip, the different method is used to search and form the shape function. Therefore, the displacement, strain, and stress of the nodes on the integral path are calculated using the Eqs. (2), (6) and (7), respectively. Combined with the SIF calculation method in Section 2.2, the Jintegration and SIF of the corresponding crack can be obtained.

q = 7.85 g/cm3 [31]. The CT sample had a size of 75  72  5 mm, as shown in Fig. 7. The fatigue pre-cracking of the specimen was performed before the test so that the total crack length was about 35 mm, the error should be kept within the range of ±0.25 mm, which can be used to study the SIF variation of the material under load. According to the reference [32], the strain gauge was used to measure the strain near the crack. Therefore, the strain gauge G1 was instrumented at the corresponding position of the crack. The orientation of the strain gauge approximately was h = 54°, r = 12.5 mm, / = 68° To study the effect of different crack configurations, two groups of samples were selected with thicknesses of 5 mm, 6 mm and 8 mm respectively. Laser cracking cracks were 30 mm and the crack width was 0.3 mm. Fatigue crack of 5 mm was generated by a fatigue load. The total crack length of the CT specimen was 35 mm, and the horizontal length from the crack tip to the load line was a = 20 mm. To reduce the environmental effect, the VIC-3D calculated the average of 5 sets of photos at a time. The camera had a focal length of 105 mm and a speckle size of 0.18 mm. The test set up as shown in Fig. 8, the fatigue and tensile loads were applied by the 25 T-MTS fatigue tester. The specimens were clamped on the fatigue tester by U-shaped clamps. The test results of the strain gauges were collected by the HBM QuantumX comprehensive data acquisition system. The two cameras of the VIC-

4. Static tensile test of CT specimen 4.1. Test set up To study the accuracy and effectiveness of the proposed hybrid DIC–EFG method, the fracture properties of the Q235qE steel compact tensile specimen (CT specimen) were studied. The mechanical properties of steel Q235qE were as follows: elastic modulus E = 210 GPa, Poisson’s ratio m = 0.33, yield strength of the materialrs = 235 MPa, tensile strengthrt = 390 MPa, density

Fig. 7. CT specimen.

6

Z. Zhu et al. / Measurement 154 (2020) 107498

Fig. 8. Test set up.

3D non-contact full-field strain measurement system was mounted at a distance of 1.2 m from the surface of the sample, the camera height was the same. The angle between the cameras was about 40°. The corresponding lighting equipment and image acquisition computer was also shown in the figure (Fig. 9). The experiment was divided into two parts, the fatigue precracking and quasi-static tensile tests. As shown in Table 1, For fatigue pre-cracking, load ratio R = 0.1 means maximum/minimum load, and the SIF at the maximum load should be less than 45% of the fracture performance (KIC) [14]. In the quasi-static tensile tests, a non-reflective black speckle pattern was applied to one side of the polished sample, and a strain gauge was attached to the other side. When the load was applied, considering the crack closure phenomenon, a preload of P0 = 100 N was applied to the specimen, and then loaded in a series of quasi-static cycles Pi, which increases step by step until the fracture of the specimen, and the load per stage was 1 KN. The difference in stress intensity factors between Pi and P0 is the result Ki. The DIC recorded the reference state of the applied preloading and the loading stated at the end of each cycle as the deformation result of the loading. 4.2. Experimental results Fatigue pre-cracking is performed according to reference [14]. Since KIC of the Q235qE material is 191.4 MPa∙m1/2 [32], when the thickness t is equal to 5 mm, according to the formula (8), the maximum load is 11.8 KN. The load amplitude is taken as 10

Fig. 9. P52 Fatigue crack growth.

KN. The fatigue pre-cracking finally causes the pre-fabricated fatigue crack length to be 5 ± 0.25 mm to stop the test. In this test, the fracture properties of samples having different thicknesses are studied. The results of fatigue pre-cracking of each group of samples are shown in Table1. From Table1, it can be seen that the fatigue pre-cracking results and the required crack calculation length a = 20 mm. The error is less than 1%, and the static tensile test can be further performed (Table 2). 4.2.1. Full-field displacement Fig. 10 shows the results of the full-field displacement test under the 10KN load of the P51 specimen. To ensure the consistency of the test results, a unified Cartesian coordinate system must be set before the calculation. The X-axis is set at the axis of symmetry. In Fig. 10 (a), under the two-way tensile load, the Xdirection displacement is shown to be symmetry related to the X-axis. The displacement reaches the maximum at the crack tip and decreases away from the axis of symmetry. In Fig. 10 (b), the displacement has obvious anti-symmetry related to the X-axis. The deformation becomes larger from left to right as it is getting closer to the load location. The displacement gradient adjacent to the crack is large. 4.2.2. Full-field strain The VIC-3D instrumentation is based on the finite element theory and aims to establish a triangular mesh model to calculate the strain. Fig. 11 shows the results of the full-field strain test under the 10KN load of the P51 specimen. The full-field strain has a symmetric distribution. And there is a concentrated three-direction strain at the crack tip. An obvious compression phenomenon is shown along the X-axis away from the crack side. The shear strain is mainly generated on both sides of the crack as well as the middle of the upper and bottom edge of the plate. The plate doesn’t have a significant deformation at the crack tip. It can be seen from the eyy cracked amplification contour that when searching with a relatively large subset, the strain around the crack and the edge of the sample could be completely missing, and the calculation results would have some effects on the strain field of the crack tip. The VIC-3D calculation indicates that when subset = 59, the fullfield average uncertainty interval reaches a minimum of only 0.01

7

Z. Zhu et al. / Measurement 154 (2020) 107498 Table 1 Experimental setup. Fatigue pre-cracking Load 10-16kN

Quasi-static tensile tests

Load ratio 0.1

Frequency 10 Hz

End condition Crack = 5 ± 0.25 mm

Preload 100 N

Load Increases per stage was 1 kN

End condition Spacimen damage

Table 2 Fatigue pre-crack test cases and results. Specimen No.

Load ratio R

Maximum load Pmax

Cycles

Fatigue crack length

Crack length a

Error

P51 P52 P61 P62 P81 P82

0.1 0.1 0.1 0.1 0.1 0.1

10 10 12 12 16 16

10,000 10,000 11,347 10,000 10,156 11,000

5.102 5.012 5.132 4.822 4.897 5.163

20.102 20.012 20.132 34.822 19.897 20.163

0.51% 0.06% 0.66% 0.89% 0.52% 0.82%

KN KN KN KN KN KN

mm mm mm mm mm mm

mm mm mm mm mm mm

Fig. 12. Validation of the strain measured by DIC at Point G1. Fig. 10. P51 full-field displacement.

pixels. To validate the test results, Fig. 12 shows the displacementload curve at the Point G1 from either strain gauge or DIC. A very good overlay is shown between the two curves. If ignoring the system error caused by the VIC-3D instrumentation, the maximum error of the measurement between the DIC and the strain gauge is 8.33%. When the load is greater than 13KN, nonlinear deformation occurs at Point G1. Thus, when subset = 59, the strain measured by the DIC is proved to be accurate by the data from the strain gauge. 4.2.3. Calculation results when using different subsets To obtain the calculated strain at the crack tip for the calculation of SIF, a smaller subset is needed for the calculation. And the

crack tip position and crack information are obtained by searching a small subset. Fig. 13 shows the DIC displacement and strain measurement results under different subsets. It can be seen that the displacement results are unchanged, but the strain results are very noisy. To study the variation of the nodes around the crack, Fig. 14 (a) shows, under different subset conditions, the displacementload curve at the crack tip G2 (12, 2) (shown in Fig. 11). The displacement increases linearly with the increase of the load. When the load is over 11KN, the sample begins to show a nonlinear deformation. Taking the test result of subset = 59 as a reference, when the subset = 39 or 19, the relative error of the displacement at Y-direction is only 0.10% or 0.17%, respectively. The displacement is shown to be the same with a different subset. Fig. 14(b) shows the strain of Point G2 with different subsets. As the subset decreases, the Y-direction strain also decreases, and even a

Fig. 11. P51 specimen full-field strain.

8

Z. Zhu et al. / Measurement 154 (2020) 107498

Fig. 13. DIC test results under different subsets.

Fig. 14. Effect of different subsets at Point G2.

negative value appears. This indicates that a small subset would lead to poor strain calculation. 5. Validation of the hybrid DIC–EFG method In this section, the meshless calculation model is developed according to the theory of Section 3 and then the full-field stress and strain are derived. The SIF calculated by the hybrid DIC-EFG method is compared with the ASTM standard. 5.1. Full-field strain based on hybrid DIC–EFG method A meshless model of CT specimens is developed based on the theory in Section 3.2, as shown in Fig. 15 (a). The VIC-3D can derive evenly distributed node information by setting the range and spacing. Step A First, DIC uses subset = 59 to calculate and uniformly distribute the field nodes according to the node spacing of

delta  delta = 2 mm  2 mm in the whole space. The transition area A´ B ´ C ´ D´ node range is x2(1.5, 15) mm, y2(11.5, 11.5) mm, and the node spacing is delta/2 = 1.0 mm. The transition region is analyzed with a small subset = 19, and the ABCD derived nodes of the inserted cracked-edge encryption region are x2(3.75, 15) mm, y2(5.75, 5.75) mm, and the node spacing is delta/4 = 0. 5 mm. Based on the EFG method, the full-field strain distribution is calculated. Fig. 15(a) shows the Y-direction strain contour of the DIC– EFG calculation on the P51 sample under 10 kN load. It can be seen from Fig. 15 (a) that the crack strike is roughly parallel to the Xaxis, so the strain at the crack tip is not symmetric about the Xaxis. The crack tip coordinate is (9.75, 0.25). Since the DIC calculation loses the high-stress gradient part of the crack tip, the strain range of the DIC test is only (0.84, 1.86)me (as shown in Fig. 11 (b)). Comparably, the DIC–EFG method can have a strain range of (1.85, 7.53)me. The hybrid DIC–EFG method can obtain the detailed information of crack tip and crack and thus can better

Fig. 15. eyy DIC–EFG calculation for P51 sample.

Z. Zhu et al. / Measurement 154 (2020) 107498

characterize the stress-strain distribution of crack tip. To validate the strain calculated by the proposed method, Fig. 15(b) shows the strain curve at Point G2 with different loads, obtained either by hybrid DIC–EFG method or DIC. Even if the scattered spots of the crack tip portion of the meshless model are composed of two parts of the DIC calculation results (subset = 59 and 19), the strain of the crack tip is still effective. The good agreement between the two curves indicates that the full-field strain calculated by the DIC–EFG method is accurate. 5.2. SIF calculation based on hybrid DIC–EFG method To calculate SIF, the integration path must be determined first. As the J-integral method has nothing to do with the integration path, the data collecting path can be far away from the crack tip high plastic area. According to Von-mises Criterion, the equivalent stress in the entire plane is calculated, and field nodes with the stress level of about 235Mpa were found to determine the plastic range. The choice of the integration path needs to include a plastic region. Finally, a 5 mm by 5 mm square with its center at the crack is selected as the integral path. (as shown in Fig. 16) Considering the crack closure phenomenon, Fig. 17 shows the SIF solution of the P51 specimen loaded and preloaded. The SIF result is the difference between the two. It can be seen from Fig. 17 that when the load is less than 12 kN, the SIF increases linearly with the increase of the load. Under 12 kN loading and pre-

Fig. 16. Von-mises stress of P51 specimen under 10kN load.

9

loading conditions, the SIF increased rapidly because of significant plastic deformation in the range of the integration path. Compared with the method that directly calculates the SIF by DIC result and Jintegral, the hybrid DIC-EFG method obtains the results around the crack, and the calculation results are more accurate. Other specimens were analyzed using the same method as the P51 specimen. Fig. 18 shows the SIF comparison between three sets of samples using either hybrid DIC–EFG method or ASTM standard solution (such as Eq. (8)). As can be seen from Fig. 18, the SIF gradually increases with the load increases, and the DIC–EFG has a similar trend with the ASTM. The DIC–EFG solution can well characterize the plastic deformation of the crack tip. Compared with the ASTM standard, the hybrid DIC–EFG method calculates the SIF based on the J-integral method. There is no need to know the load and boundary conditions, nor the crack configuration or the geometry of the test sample, only the displacement field and strain field adjacent to the crack need to be measured.

5.3. Sensitivity analysis For DIC, the calculation results of SIFs under different subset conditions are analyzed. The DIC-EFG calculation model consists of two parts, S0 and the crack tip field node Sc (as shown in Fig. 5). The subset calculated in S0 satisfies the minimum error of image matching, which is given by the DIC instrument. The subset of the crack tip Sc needs to satisfy the search for the crack boundary. Therefore, the subset of the crack tip need more analysis. Fig. 19 shows the variation of SIF with load in 6 different subsets of 0–59, while subset = 0 means that there is no Sc part in the calculation model. It can be seen that, compared with ASTM standard, the results of subset = 0, 49 and 59 have large errors. When subset = 49, the maximum error between the SIF calculated by the mixed DIC-EFG method and the ASTM standard solution reached 26.8%. It is because large subset or without crack-tip field node will missing the strain around the crack. When calculating J-integral, the nodes cannot form a closed fence, which results in a large error SIF value. When the subset is less than 39, the results of the three tests are basically the same, and the error between the SIF calculated by the mixed DIC-EFG method and the ASTM standard solutions are less than 11.9%. Therefore, when subset = 19, the hybrid DIC-EFG method has better stability. Considering the influence of node density on the calculation model, the sensitivity of the SIFs results under the condition of delta = 1–4 is studied. It can be seen from Fig. 20 that when delta = 3,4, most of the SIFs results are larger than the ASTM standard with a maximum error of 32.3%. When delta = 1,2, the SIFs results tend to be stable, and the maximum error compared with

Fig. 17. SIF Solution of P51 Specimen.

10

Z. Zhu et al. / Measurement 154 (2020) 107498

Fig. 18. SIF calculation results.

Fig. 19. SIFs under different subset conditions.

the ASTM standard is 12.9%. so the final model node spacing of delta = 2 mm is able to meet the calculation accuracy requirements. For the EFG method, the accuracy of meshless calculation depends on the nodes in the influence domain, so choosing an appropriate influence domain can improve the calculation accuracy. For a calculation point, the size of the support domain di should be larger than the delta between nodes. Liu [33] shows that in general, taking (2.0 ~ 3.0)  delta for most problems can achieve better analysis results. Aiming at the research content of this paper, and taking the samples of P51 as examples, the influence of the size of the influence region in the range of (1.5 ~ 4.0)  delta on the SIF calculation results was studied. From the sensitivity analysis of di parameters, it can be seen from Fig. 21 that with the increase of the load, the variation of SIFs under different di is basically the same. Taking the calculation result of 10 kN as an example, it can be seen that the error of SIFs between different di values and the ASTM standard within ± 5%. The value of di parameter in this article is 3, which can be con-

Fig. 20. SIF solution under different delta.

Fig. 21. di parameter sensitivity analysis P51 specimen.

cluded that the hybrid DIC-EFG method proposed in this paper has good stability with the support domain di.

Z. Zhu et al. / Measurement 154 (2020) 107498

6. Conclusion This paper proposes a hybrid DIC–EFG method for extracting the SIF of the surface crack. The DIC is used to obtain the field node information and the full-field displacement. The EFG method can directly calculate the full-field stress and strain, and then the Jintegral method is used to calculate the SIF of the fatigue crack. Using a VIC-3D, a quasi-static cyclic static tensile test is conducted on a compact tensile specimen (CT) with Q235 material to validate the proposed method. The conclusions are as follows: (1) VIC-3D is used to retrieve the displacement field with good stability, and it is stable with the different selected subset, and the calculated strain is greatly affected by the size of the subset. (2) Based on DIC, the method of deriving nodes to develop a meshless model does not need to be remodeled and meshed for the crack tip. The modeling is simple for structures having irregular geometry and cracks. The full-fie‘ld measured displacement from DIC is introduced into the meshless model as a boundary condition. This method is easy to implement without knowing the load and constraints. The hybrid DIC–EFG method can obtain the accurate full-field strain including the crack tip, which is a great challenge for the DIC. (3) Through the study of the quasi-static static tensile test on the CT specimens and parameter sensitive analysis, the hybrid DIC–EFG method is validated to be accurate and stable. The results show that the method can effectively characterize the actual stress-strain state of the crack tip and calculate the SIF. The proposed method can be used for fracture and fatigue properties of cracked structures in the future.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments Financial support for this study was provided by the National Natural Science Foundation of China (Grant No.51678576), the National Key R&D Program of China (Grant No. 2017YFB1201204) and the Fundamental Research Funds for the Central Universities of Central South University (Grant No. 1053320182867). References [1] H. Li, B. Zhao, H. Zhu, Numerical Simulation of Fatigue Performance of Diaphragm of Large-Span Bridge Orthotropic Deck, Complexity (2018). [2] Y. Chen, B. Zhang, N. Zhang, M. Zheng, A condensation method for the dynamic analysis of vertical vehicle–track interaction considering vehicle flexibility, J. Vibr. Acoust. 137 (2015) 041010. [3] G. Yi, T. Yu, T.Q. Bui, C. Ma, S. Hirose, SIFs evaluation of sharp V-notched fracture by XFEM and strain energy approach, Theor. Appl. Fract. Mech. 89 (2017) 35–44. [4] C. Hwu, H.Y. Huang, Investigation of the stress intensity factors for interface corners, Eng. Fract. Mech. 93 (2012) 204–224.

11

[5] C. Roux-Langlois, A. Gravouil, M.C. Baietto, J. Réthoré, F. Mathieu, F. Hild, et al., DIC identification and X-FEM simulation of fatigue crack growth based on the Williams’ series, Int. J. Solids Struct. 53 (2015) 38–47. [6] Z. Fu, B. Ji, C. Zhang, D. Li, Experimental study on the fatigue performance of roof and U-rib welds of orthotropic steel bridge decks, KSCE J. Civ. Eng. 22 (2018) 270–278. [7] S. Roux, J. Réthoré, F. Hild, Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2D and 3D cracks, J. Phys. D Appl. Phys. 42 (2009) 214004. [8] Q. Tang, J. Hu, T. Yu, Electromagnetic evaluation of brick specimens using synthetic aperture radar imaging, NDTE Int. 104 (2019) 98–107. [9] C. Dong, X.W. Ye, T. Liu, Non-contact structural vibration monitoring under varying environmental conditions, Vibroengineering Procedia 5 (2015) 217– 222. [10] B.V. Farahani, P.J. Tavares, J. Belinha, P. Moreira, A fracture mechanics study of a compact tension specimen: digital image correlation, finite element and meshless methods, Procedia Struct. Integrity 5 (2017) 920–927. [11] Hadadian A, Afaghi Khatibi A, Nikkhah M, Mohammadi F. A methodology for evaluating mixed-mode (I+ II) SIF, combining experimental data acquired by ESPI with Westergaard type solution. World Congress on Engineering 2010: International Association of Engineers; 2010. p. 1251-1256. [12] X. Dai, Q. Pu, L. Wang, H. Yun, Y. Wang, Measurement on fracture process and prediction of the load capacity of steel fiber reinforced concrete by electronic speckle pattern interferometry, Compos. B 42 (2011) 1181–1188. [13] Dong C, Celik O, Catbas F, OBrien E, Taylor S. A Robust Vision-based Method for Displacement Measurement under Adverse Environmental Factors using Spatio-Temporal Context Learning and Taylor Approximation. 2019. [14] S.M. Barhli, M. Mostafavi, A.F. Cinar, D. Hollis, T.J. Marrow, J-integral calculation by finite element processing of measured full-field surface displacements, Exp. Mech. 57 (2017) 997–1009. [15] R. Zhang, L. He, Measurement of mixed-mode stress intensity factors using digital image correlation method, Opt. Lasers Eng. 50 (2012) 1001–1007. [16] F. Zhong, C. Quan, Efficient digital image correlation using gradient orientation, Optics and Lasers Technology. 106 (2018) 417–426. [17] S.G. Shah, J.M.C. Kishen, Determination of fracture parameters of concrete interfaces using, DIC (2010). [18] T. Fawad, S.M. Zeeshan, N. Nausheen, Practical Application of DIC in Fatigue and Fracture Toughness, Testing. (2012). [19] Y. Chen, P. Logan, P. Avitabile, Non-model based expansion from limited points to an augmented set of points using Chebyshev polynomials, Exp. Tech. 1–23 (2019). [20] Chen Y, Joffre D, Avitabile P. Underwater dynamic response at limited points expanded to full-field strain response. 2018;140:051016. [21] C.K. Desai, S. Basu, V. Parameswaran, Determination of complex stress intensity factor for a crack in a bimaterial interface using digital image correlation, Optics Lasers in Eng. 50 (2012) 1423–1430. [22] M.S. Reddy, K. Ramesh, A. Thiyagarajan, Evaluation of mode-I SIF, T-stress and J-integral using displacement data from digital image correlation – Revisited, Theor. Appl. Fract. Mech. 96 (2018) 146–159. [23] S.R. Mcneill, W.H. Peters, M.A. Sutton, Estimation of stress intensity factor by digital image correlation, Eng. Fract. Mech. 28 (1987) 101–112. [24] S. Roux, F. Hild, Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches, Int. J. Fract. 140 (2006) 141–157. [25] K. Bharadwaj, A. Sheidaei, A. Afshar, J. Baqersad, Full-field strain prediction using mode shapes measured with digital image correlation, Measurement (2019). [26] G.L. Gonzáles, J.A. González, J.T. Castro, J.L. Freire, A J-integral approach using digital image correlation for evaluating stress intensity factors in fatigue cracks with closure effects, Theor. Appl. Fract. Mech. 90 (2017) 14–21. [27] G. Gonzáles, M. Meggiolaro, Strain field measurements around notches using SIFT features and meshless methods, Appl. Opt. 54 (2015) 4520–4528. [28] M. Hajiazizi, P. Bastan, The elastoplastic analysis of a tunnel using the EFG method: A comparison of the EFGM with FEM and FDM, Appl. Math. Comput. 234 (2014) 82–113. [29] J.R. Rice, A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks, J. Appl. Mech. 35 (1968) 379–386. [30] G.L.G. Gonzáles, J.A.O. González, J.T.P. Castro, A J-integral approach using digital image correlation for evaluating stress intensity factors in fatigue cracks with closure effects, Theor. Appl. Fract. Mech. 90 (2017) 14–21. [31] S. Guangyong, L. Xinglong, Z. Gang, On fracture characteristics of adhesive joints with dissimilar materials – An experimental study using digital image correlation (DIC) technique, Compos. Struct. 201 (2018) 1056–1075. [32] Z. Zhangyan, L. Yunbing, S. Guozheng, Experiment Measuring Fracture Toughness of Q235 Steel by J Integral, J. Wuhan Univ. Technol. 24 (2002) 111–112. [33] G.R. Liu, Y.T. Gu. An introduction to meshfree methods and their programming. 2005.