Combing DGS and finite element for stress analysis using inverse boundary method Rui Zhang,1,* Ran Guo,1,2 and Heming Cheng1,3 1
Engineering Mechanics Department of Kunming University of Science and Technology, Kunming 650500, China 2
e-mail: [email protected]
3
e-mail: [email protected]
*Corresponding author: [email protected] Received 8 August 2014; revised 22 September 2014; accepted 3 October 2014; posted 11 November 2014 (Doc. ID 220642); published 10 December 2014
Digital gradient sensing (DGS) data combined with the finite-element method is proposed for stress solutions over the stress concentration area. Boundary conditions for a local finite element model, that is, the nodal force along the boundaries, are inversely determined from experimental values obtained by the DGS method. The DGS method measures the Cartesian stress gradient components directly. The sum in Cartesian stresses at all interesting points on the surface is obtained from the stress gradient using the linear least squares method. Thus, the sum stresses are used to compute the unknown boundary conditions for the local model. After boundary conditions are computed, the individual stress components are calculated by the direct finite element method. The effectiveness is demonstrated by applying the proposed method to a three-point bending specimen under the compression problem. Results show that the boundary conditions of the local finite element model can be determined from the DGS data and then the individual stresses can be obtained by the proposed method. © 2014 Optical Society of America OCIS codes: (200.0200) Optics in computing; (200.4560) Optical data processing; (350.4800) Optical standards and testing. http://dx.doi.org/10.1364/AO.53.008350
1. Introduction
Digital image correlation (DIC) is a practical and effective optical technique for surface deformation measurement. By processing the digital images of a test object surface recorded by an imaging device before and after deformation, DIC directly provides full-field displacements of the object surface with subpixel accuracy, and strains can be calculated directly or by further manipulations of the displacement fields. DIC has been widely used for full-field deformation measurement in the field of experimental mechanics due to its advantage of full-field measurement, easy specimen preparation, simple optical arrangement, and experimental setup.
1559-128X/14/358350-08$15.00/0 © 2014 Optical Society of America 8350
APPLIED OPTICS / Vol. 53, No. 35 / 10 December 2014
DIC methods give objective displacement information over the whole field. It should be noted that strains obtained from the N-R method are not accurate. So a lot of effort has been made to improve the accuracy of the strain. The gradients of displacement (strain fields) can be calculated by a direct finite difference approximation, which will amplify the noise and give a poor result in experimental cases; this is especially true if the strains around the load application point are of interest. For this reason, many differentiation algorithms have been developed to calculate strain fields from noisy displacement fields. Some finite element techniques for smoothing and differentiating the experimental displacement fields are proposed [1–4]. The measured displacements are used as the displacements at the nodes of elements of a finite element model. Then, the strains inside the elements are calculated using a displacement-strain
relationship. The disadvantage of the techniques utilizing the finite element method (FEM) is the preparation of a finite element model. Another approach to the calculation of the strains from the measured displacements is the least squares method. In this case, a polynomial is fitted to displacement values in a small local region by a least squares method. The strains in the region are then obtained from the coefficients of the approximated polynomial. As a result, the noise of the displacement field could be reduced and the strains could be exacted. However, the high-frequency part of the displacement is also eliminated and this is particularly harmful in the case of a high displacement gradient. The measured strains at the boundary also frequently contain relatively large errors compared with those inside the boundary. So it is meaningful to have a method that offers all the advantages of DIC yet is capable of directly measuring stress gradients in the whole field near stress concentrators. Accordingly, as an extension of the DIC-based method, authors proposed a full-field optical technique called Digital Gradient Sensing (DGS) [5,6] which can quantify the elastooptic effect using the DIC method for mechanical characterization of optically transparent planar solids. The method can link angular deflections of light rays quantified using DIC to two orthogonal stress gradients under plane stress conditions. So, the DGS measures Cartesian stress gradient components directly; subsequently, Chandru Periasamy [5] calculated the sum in Cartesian stresses or radial stress in the special case. Significantly, sum stress can be obtained by numerical integration from the stress gradient. Using the integral method has the following advantages compared with the differential method. Differential has an amplification effect on the high frequency signal. For example, it usually needs filtering and smoothing in the process of image processing while integral suppresses the high frequency part, and has no special requirements on the smooth processing. However, individual stress components that directly relate to the material mechanical properties are preferred, and a stress separation method is essential. This is a similar issue as that for photoelasticity, which measures the principal stress difference and the principal direction, and thus the stress components themselves cannot be obtained directly. Several methods have been proposed for stress separation in photoelasticity such as the shear difference method [7,8]; numerical analysis of the compatibility equation, oblique incidence and hybrid methods of photoelasticity, and another experimental method [9–12]. However, all of the above methods for stress separation in photoelasticity are not suitable for that in the DGS method because shear stress is not available when the shear difference method implants in the DGS method. Hybrid experimental methods require a significant modification in the experimental setup, which is a cause of inconvenience, and the other two methods are also impracticable due to their
difficult measurement theory. On the other hand, several hybrid methods of photoelasticity and numerical analysis [13–15], and inverse analysis methods have also been proposed for the stress separation. D. G. Berghaus [13] provides a hybrid method with photoelasticity and finite-element combined, from which displacements along free surfaces and axes of symmetry are obtained from photoelasticity and then perform the FEM to get stress components. S. Yoneyama [15] also provides a hybrid method with photoelasticity and finite-element combined, from which tractions along boundaries are obtained from both principal stress difference and principal direction, and then perform the FEM to get the stress components. In this study, an inverse method is proposed for stress separation in the DGS. First, the stress sum of each point is estimated from the stress sum gradient using the linear least squares method. Then, the interesting area is subdivided into quadrilateral elements, nodal forces along boundary are obtained by the stress sum in the DGS and individual stress components are obtained by the FEM directly. Finally, the inverse method is applied to a three-point bending specimen under compression test. Stress components around the crack tip are obtained; stress components at each point are compared with the FEM direct analysis for this problem; and the results are summarized and conclusions are drawn. 2. Theoretical Background A. Experimental Procedure
Figure 1 schematically illustrates the typical experimental setup for the DGS. The experiment contains a speckle target, a transparent test specimen, and a digital camera. The target plane surface must have a random speckle pattern. The transparent specimen to be tested is placed in front of and parallel to the target plane at a distance, Δ. The target plane is uniformly illuminated using two white light sources. The camera is placed in front of the test specimen at a long distance, L (L ≫ Δ), with its optical axis normal to the specimen surface, imaging the target plane surface through the specimen in the region of interest. The image captured under no-load is
Fig. 1. Experimental setup for DGS. 10 December 2014 / Vol. 53, No. 35 / APPLIED OPTICS
8351
referred to as a reference image or undeformed image. The image captured under loading is referred to as a deformed image. The speckles on the target plane are compared to the undeformed image due to the local change of the thickness and refractive index at a load. So displacements recorded using DIC can be related to angular deflection fields. These angular deflections can further be related to the inplane stress gradient, as discussed in the following section. The measurement principle of the DGS as shown in Fig. 2 is as follows. The thickness of a transparent specimen is B and refractive index is n0 in a no-load state. Cartesian systems, x; y; z and x0 ; y0 ; z0 , are chosen on the specimen and target plane, respectively, so that the z axis is perpendicular to the specimen (or target plane). When the target plane surface is imaged through a transparent specimen, a point, P, on the target plane corresponding to O on the specimen is imaged in the reference state. When imposing a load, the thickness and refractive index of a specimen are all changed due to local stress; thus, a neighbor point, Q, corresponding to O is imaged in the deformed state. That is, light rays, OP, in a reference state corresponding to OQ in a deformed state. The spatial vector, PQ, with known Δ can be related to angular deflections of light rays. On the assumption of isotropic, linearly elastic, and plane stress conditions, the governing equations of the DGS method are expressed as [5] ∂σ xx σ yy δx ; ≈ α CB ∂x Δ δy ∂σ xx σ yy φy ≈ β CB : ∂y Δ φx
(1)
From above, δx and δy can be derived from DIC, Δ is measured in advance, and C and B as material property are also known. So, angular deflections as well as stress gradients are obtained. Stress sums σ xx σ yy , determined from measured stress gradients, are shown in Section 4.
B. Inverse Boundary Method Using Cartesian Stress Components
The interesting region for stress analysis is subdivided into quadrilateral elements as shown in Fig. 3. For a finite element model, if the boundary conditions are applied appropriately, then stress components at every interesting point could be obtained accurately. However, for the local finite element model matching the analysis area, the boundary conditions are usually unknown. The unknown boundary nodal forces around the analysis region are inversely calculated using the measured sums of the Cartesian stress components at the interior interesting points. Then, individual stress components at interior points are determined by finite element direct analysis using the computed boundary conditions. The basic theory of the hybrid method is shown in Fig. 3. The displacements of two orthogonal directions at node, c1 , and the displacements of the y direction at node, c2 , are fixed, so rigid body motion is not allowed. Apply a unit force on an element node along the boundary with its direction in accordance with the x or y axis. Do the finite element procedure and then the stress components at any interesting point can be calculated at this applied unit force. The procedure is repeated by changing the applied node and direction. That is, for the applied load of Pj 1Nj 1…n, stress components of σ 0x ji ; σ 0y ji i 1…m at any interesting data point such as Qj xi; yi are obtained, where j is the unit force index, i is data point index, n is the total number of the unknown nodal forces to be determined, and m is the total number of interesting point. Thus, the stress components of σ x ji and σ y ji at any point xi; yi under real nodal force can be expressed as follows: σ x i σ y i
n X j1 n X j1
σ 0x ji F j i 1…m; σ 0y ji F j i 1…m;
(2)
where F j is the real nodal force along the boundary.
Fig. 2. Schematic of the working principle of the DGS. 8352
APPLIED OPTICS / Vol. 53, No. 35 / 10 December 2014
Fig. 3. Finite element model with unit force applied on the boundary.
So the sum of the stress components (stress sum) obtained by the DGS method can be related to unknown nodal force as follows: σ x σ y i
n X σ 0x ji σ 0y ji F j i 1…m;
(3)
fatigue loading in pure mode I at 1 Hz, average load of 1 kN, and alternating load of 0.3 kN, resulting in a total crack length of 12.2 mm. Then, black paint was sprayed on the target plane surface in order to apply a high contrast feature on the specimen surface.
j1
B. Sum Stress Gradient Measurement
where σ x σ y i is the stress sum obtained by the DGS at the point xi; yi under an unknown applied load of F j. Eq. (3) can be written in matrix form as T bF;
(4)
where 2 6 T6 4 2 6 b6 4
σ x σ y 1 .. .
3 7 7 5
σ x σ y m σ 0x 11 σ 0y 11 .. .
σ 0x 1m σ 0y 1m 3 F1 6 . 7 7 F6 4 .. 5: 2
..
σ 0x n1 σ 0y n1 .. .
.
3 7 7 5
σ 0x nm σ 0y nm
Fn The solution of Eq. (4) in the least squares sense is F bT b−1 bT T:
(5)
Equilibrium equations relate the total applied loads in the horizontal and/or vertical directions to the loads produced by stresses at experimental data points. Using the above equations, the applied loads (nodal forces) along the boundary are calculated. Stress components at every experimental data point are then computed under the calculated nodal forces.
The PMMA specimen is placed in front of a target plane painted with a black and white speckle pattern with (20 mm) distances. A CCD is placed at a large distance (L ≈ 680 mm) from the specimen is combined with a computer for capturing images. At a small load of several Newtons, the image is recorded using a camera (2048 × 2448 pixels) as a reference image. Subsequently, speckle images at different load levels were captured as deformed images. As fields around the crack tip suffer deformation, light rays passing through the specimen were deflected by nonuniform stress, thus the speckle image was deformed relative to the reference state. Displacements were analyzed by the DGS, so full-field displacements and, subsequently, angular deflections (stress gradients) were extracted. During the DGS analysis of these images, an area of 7.69 mm × 5.77 mm around the crack tip is chosen to be the region of interest, as shown in Fig. 4. In the image, 1 mm corresponds to 26 pixels. Angular deflections were calculated with subset size of 29 × 29 pixels and grid step of 5 pixels, which yielded a total number of 1271 with an array of 31 × 41 data points. When calculating stress components we consider rigid body motions and imposing boundary conditions of the problem to quantify the analysis of the contour levels. That is, in the current problem, the boundary conditions, such as asymmetric stress gradients (φy ) in the y direction about the x axis, symmetric stress gradients (φx ) in the x direction about the x axis, and vanishing stress gradients far away from the loading point are all utilized. Thus, Fig. 5 shows the resulting stress gradient of ∂σ xx σ yy ∕∂x and ∂σ xx σ yy ∕∂y with the rigid body translation deleted in the end for the representative load of 600 N in a square region around the crack tip, and the crack tip position is (30.5 mm, 34.5 mm).
3. Experimental Procedure A.
Specimen and Equipment
To verify the proposed hybrid method, a stress concentration problem was studied. An edge-notched three-point bending specimen (220 mm × 60 mm × 9 mm) was manufactured from polymethyl methacrylate (PMMA). The elastic modulus is 2.8 GPa and Poisson’s ratio is 0.37, which are measured by strain gauges; and the stress optical constant is −1.0 × 10−10 m2 ∕N, which calibrated in the Appendix described in the literature [16]. In order to simulate that the crack exists in the actual components, the edge-notched PMMA specimen was subjected to
Fig. 4. Experimental analysis area. 10 December 2014 / Vol. 53, No. 35 / APPLIED OPTICS
8353
Fig. 5. Contours of the sum stress gradient for the x direction (left) and the y direction (right) for 600 N.
C.
2
Sum Stress Calculated by Linear Least Squares
Numerical integration converts gradient data into sum stress by using least squares fitting, pi;j1 pi;j T i;j1 − T i;j ; 2 xi;j1 − xi;j i 1; …; M;
J 1; …; N − 1;
qi1;j qi;j T i1;j − T i;j ; 2 yi1;j − yi;j i 1; …; M − 1;
J 1; …; N;
where x, y, are the Cartesian coordinates, p and q denote the stress gradients in the x and y directions, respectively, i.e., p ∂σ xx σ yy ∕∂x and q ∂σ xx σ yy ∕∂y. T is the stress sum, i.e., T σ xx σ yy . M and N are the total rows and total columns about the calculated area. According to Eq. (6), the stress sum information can be acquired as 1 T i;j1 − T i;j pi;j1 pi;j xi;j1 − yi;j ; 2 i 1; …; M; J 1; …; N − 1; 1 T i1;j − T i;j qi1;j qi;j xi1;j − yi;j ; 2 i 1; …; M − 1; J 1; …; N:
(7)
Furthermore, Eq. (7) can be rewritten in terms of matrices as
6 6 6 T 6 6 6 4 2
where 8354
APPLIED OPTICS / Vol. 53, No. 35 / 10 December 2014
0 −1 1 3 T 1;1 7 T 2;1 7 7 7 ; .. 7 . 7 5
T M;N
;
M−1NMN−1×MN
MN×1
3 p1;2 p1;1 x1;2 −x1;1 6 7 6 p1;3 p1;2 x1;3 −x1;2 7 6 7 6 7 .. 6 7 6 . 7 6 7 6 p 6 M;N pM;N−1 xM;N −xM;N−1 7 7 G 6 7 q2;1 q1;1 y2;1 −y1;1 6 7 6 7 6 7 q y −y q 3;1 2;1 3;1 2;1 6 7 6 7 6 .. 7 6 . 5 4 qM;N qM−1;N yM;N −yM−1;N
:
M−1NMN−1×1
Least squares estimation can be carried out to provide the stress sum as an optimized solution. The counters of the stress sums thus obtained are shown in Fig. 6. D.
AT G;
3
6 0 −1 0 0 1 0 0 7 7 6 6 . . . . . . . . .7 6 . . . . . . . . .7 6 . . . . . . . . .7 7 6 6 0 0 −1 0 0 1 7 7 6 A 6 7 6 −1 1 0 0 7 7 6 6 0 −1 1 0 0 7 7 6 6 . . . . . . . . .7 6 . . . . . . . . .7 4 . . . . . . . . .5 2
(6)
−1 0 0 1 0 0
Stress Components Obtained by Inverse Boundary
The stress separation is performed in the 8 mm × 6 mm region around the crack tip. This area is divided into a finite element mesh which is a little larger than the experimental areas in order to
Fig. 8. Nodal force along a boundary. Fig. 6. Estimated stress sum obtained by the linear least squares method for 600 N.
include all experimental data points. Using the procedures as Section 2 described to calculate the stress for each experimental data point x; y, it must be found to which element it belongs and mapped from global coordinates x; y to local coordinator ε; η. Numerical interpolation is adopted here to obtain the discrete expressions of the mapping [17], considering the fact that it is generally not possible to determine the explicit expressions of the mapping. Figure 7 shows the finite element model of the 8 mm × 6 mm region used for the proposed method. In this model, 4-noded isoparametric elements are used. The numbers of the elements and the nodes are 300 and 336, respectively. In order to obtain the stresses under the unit force at a point on the boundary, the displacements at some nodes must be fixed to prevent the rigid body motion. Displacements of the y direction at node c1 and c2 are fixed; displacements of the x direction at node c3 are fixed. So, rigid body motion is not allowed, and other nodal force on the boundary needs to be calculated. The total number of nodes on the boundary is 70, so the nodal force to be determined is 137 because the three displacement components at the points, c1 , c2 , and c3 are fixed. The nodal forces at 137 points are simultaneously determined as described in Section 2. The nodal forces at 19 points along boundary, c4 ; c5 obtained from the estimated stress sum, are shown in Fig. 8. For comparison, the nodal force obtained from the FEM direct analysis are also shown. As
Fig. 7. Finite element model for analysis area.
shown in this figure, the significant difference between the nodal forces obtained using the proposed method and direct analysis is not observed. Stress components could then be obtained by finite element analysis directly. Obtaining the nodal force and acquisition of stress components are performed using Fortran language programs made by the authors. The numbers of finally separated stress components obtained from the experimental stress sums are shown in Fig. 9. An FEM analysis of σ x and σ y is also shown in Fig. 10. By comparison, we can see that the calculated stress component distribution by experimental stress sum corresponds well to theoretical calculations. E. Discussion of the K -Dominance Characterization About the Selected Area
The optical measurements are usually used in evaluating crack-tip parameters. Crack-tip parameters such as stress intensity factor (SIF) can be determined from full-field in-plane optical data based on the generalized polynomial functions proposed by Williams. A K-dominant field can be defined as the one in which the contribution from the higherorder terms is negligible when compared to the first term. If there is good agreement between experimental and theoretical values when only the first term in the expansion field is used, it suggests that the use of the K-dominant field itself may be sufficient in obtaining an accurate parameter in the specimen’s geometry. However, due to the finite size of the specimens, higher-order nonsingular deformations generally influence the crack-tip fields. Hence, the assumptions of K-dominance are relaxed and higher order terms are included. According to the literature [18,16], there exists a 3D zone surrounding the crack tip and a possible lack of K-dominance. The maximum extent of the 3D region corresponds to approximately half of the specimen thickness. So, the selected area of 7.69 mm × 5.77 mm around the crack tip might lack K-dominance. 3D effects generally prevent data from within r∕B 0.5 (r is the distance from the crack tip and B is the specimen thickness) from being analyzable using 2D 10 December 2014 / Vol. 53, No. 35 / APPLIED OPTICS
8355
Fig. 9. Stress components of σ x and σ y for 600 N determined by the proposed inverse method.
Fig. 10. Stress components of σ x and σ y estimate by finite element direct analysis.
descriptions. For mode I cracks as the mode in this paper, however, both optical and finite-element [19] results have revealed that there exist sectors in the crack-tip vicinity where 3D effects are minimal and the 2D singular term is dominant. The data in these regions, although well within the r∕B 0.5 limit, can be well described by the K-dominant field. So mode I fracture is dependent on the critical values of K I . 4. Conclusions
The inverse method theory for DGS stress analysis is developed using the sum of the Cartesian stresses at interesting data points, which can be calculated from the DGS method when aided by analytical expressions. The DGS measurements of the stress gradient in Cartesian stresses at all interesting points are used to compute the stress sum of each point, and stress sums are used for the inverse boundary analysis. Then the individual stress components are calculated by the FEM directly using the computed boundary conditions. The method has also been implemented to study a stress concentration problem about the crack tip. The results indicate that the actual boundary nodal forces can be predicted with a reasonable accuracy, and stress components can also be obtained reasonably. The work is supported by the National Basic Research Program of China (under Grant 14118683); 8356
APPLIED OPTICS / Vol. 53, No. 35 / 10 December 2014
China Postdoctoral Science Foundation (under Grant 2014M552548XB); and the National Natural Science Foundation of China (under Grant 11402103). References 1. L. B. Meng, G. C. Jin, and X. F. Yao, “Application of iteration and finite element smoothing technique for displacement and strain measurement of digital speckle correlation,” Opt. Lasers Eng. 45, 57–63 (2007). 2. M. A. Sutton, J. L. Turner, H. A. Bruck, and T. A. Chae, “Fullfield representation of discretely sampled surface deformation for displacement and strain analysis,” Exp. Mech. 31, 168–177 (1991). 3. D. J. Segalman, D. B. Woyak, and R. E. Rowlands, “Smooth splike-like finite-element differentiation of fullfield experimental data over arbitrary geometry,” Exp. Mech. 19, 429–437 (1979). 4. B. Pan, H. M. Xie, Z. Q. Guo, and T. Hua, “Full-field strain measurement using a two- dimensional Savitzky–Golay digital differentiator in digital image correlation,” Opt. Eng. 46, 033601 (2007). 5. C. Periasamy and H. V. Tippur, “A full-field digital gradient sensing method for evaluating stress gradients in transparent solids,” Appl. Opt. 51, 2088–2097 (2012). 6. C. Periasamy and H. V. Tippur, “Measurement of crack-tip and punch-tip transient deformations and stress intensity factors using Digital Gradient Sensing technique,” Eng. Fract. Mech. 98, 185–199 (2013). 7. S. J. Haake and E. A. Patterson, “The determination of principal stresses from photoelastic data,” Strain 28, 153–158 (1992). 8. M. Solaguren-Beascoa Fernández, J. M. Alegre Calderon, P. M. Bravo, and I. I. Cuesta Segura, “Stress-separation techniques in photoelasticity: a review,” J. Strain Anal. Eng. Des. 45, 1–17, (2010).
9. R. J. Greene and E. A. Patterson, “An integrated approach to the separation of principal surface stresses using combined thermophoto-elasticity,” Exp. Mech. 46, 19–29 (2006). 10. T. Sakagami, S. Kubo, and Y. Fujinami, “Full-field stress separation using thermoelasticity and photoelasticity and its application to fracture mechanics,” JSME Int. J. Ser. A 47(3), 298–304 (2004). 11. S. Yoneyama, Y. Morimoto, and M. Kawamura, “Twodimensional stress separation using phase-stepping interferometric photoelasticity,” Meas. Sci. Technol. 16, 1329–1334 (2005). 12. J. Lim and K. Ravi-Chandar, “Dynamic measurement of two dimensional stress components in birefringent materials,” Exp. Mech. 49, 403–416 (2009). 13. D. G. Berghaus, “Combining photoelasticity and finite element methods for stress analysis using least squares,” Exp. Mech. 31, 36–41 (1991). 14. D. Chen, A. A. Becker, I. A. Jones, T. H. Hyde, and P. Wang, “Development of new inverse boundary element techniques
15. 16. 17. 18.
19.
in photoelasticity,” J. Strain Anal. Eng. Des. 36, 253–264 (2001). S. Yoneyama, S. Arikawa, and Y. Kobayashi, “Linear and nonlinear algorithms for stress separation in photoelasticity,” Exp. Mech. 52, 529–538 (2012). R. Zhang, R. Guo, and S. Wang, “Mixed mode fracture study of PMMA using digital gradient sensing method,” Eng. Fract. Mech. 119, 164–172 (2014). Y. P. Yu, Y. Xiang, Q. S. Liu, and A. Zhong, “Research on numerical inverse isoparametric mapping interpolation and its application,” Rock Soil Mech. 22, 226–228 (2001) (in Chinese). S. Ramaswamy, H. V. Tippur, and L. Xu, “Mixed-mode cracktip deformation studied using a modified flexural specimen and coherent gradient sensing,” Exp. Mech. 33, 218–227 (1993). H. V. Tippur, S. Krishnaswamy, and A. J. Rosakis, “Optical mapping of crack tip deformations using the methods of transmission and reflection coherent gradient sensing: a study of crack tip K-dominance,” Int. J. Fract. 52, 91–117 (1991).
10 December 2014 / Vol. 53, No. 35 / APPLIED OPTICS
8357
Engineering Mechanics Department of Kunming University of Science and Technology, Kunming 650500, China 2
e-mail: [email protected]
3
e-mail: [email protected]
*Corresponding author: [email protected] Received 8 August 2014; revised 22 September 2014; accepted 3 October 2014; posted 11 November 2014 (Doc. ID 220642); published 10 December 2014
Digital gradient sensing (DGS) data combined with the finite-element method is proposed for stress solutions over the stress concentration area. Boundary conditions for a local finite element model, that is, the nodal force along the boundaries, are inversely determined from experimental values obtained by the DGS method. The DGS method measures the Cartesian stress gradient components directly. The sum in Cartesian stresses at all interesting points on the surface is obtained from the stress gradient using the linear least squares method. Thus, the sum stresses are used to compute the unknown boundary conditions for the local model. After boundary conditions are computed, the individual stress components are calculated by the direct finite element method. The effectiveness is demonstrated by applying the proposed method to a three-point bending specimen under the compression problem. Results show that the boundary conditions of the local finite element model can be determined from the DGS data and then the individual stresses can be obtained by the proposed method. © 2014 Optical Society of America OCIS codes: (200.0200) Optics in computing; (200.4560) Optical data processing; (350.4800) Optical standards and testing. http://dx.doi.org/10.1364/AO.53.008350
1. Introduction
Digital image correlation (DIC) is a practical and effective optical technique for surface deformation measurement. By processing the digital images of a test object surface recorded by an imaging device before and after deformation, DIC directly provides full-field displacements of the object surface with subpixel accuracy, and strains can be calculated directly or by further manipulations of the displacement fields. DIC has been widely used for full-field deformation measurement in the field of experimental mechanics due to its advantage of full-field measurement, easy specimen preparation, simple optical arrangement, and experimental setup.
1559-128X/14/358350-08$15.00/0 © 2014 Optical Society of America 8350
APPLIED OPTICS / Vol. 53, No. 35 / 10 December 2014
DIC methods give objective displacement information over the whole field. It should be noted that strains obtained from the N-R method are not accurate. So a lot of effort has been made to improve the accuracy of the strain. The gradients of displacement (strain fields) can be calculated by a direct finite difference approximation, which will amplify the noise and give a poor result in experimental cases; this is especially true if the strains around the load application point are of interest. For this reason, many differentiation algorithms have been developed to calculate strain fields from noisy displacement fields. Some finite element techniques for smoothing and differentiating the experimental displacement fields are proposed [1–4]. The measured displacements are used as the displacements at the nodes of elements of a finite element model. Then, the strains inside the elements are calculated using a displacement-strain
relationship. The disadvantage of the techniques utilizing the finite element method (FEM) is the preparation of a finite element model. Another approach to the calculation of the strains from the measured displacements is the least squares method. In this case, a polynomial is fitted to displacement values in a small local region by a least squares method. The strains in the region are then obtained from the coefficients of the approximated polynomial. As a result, the noise of the displacement field could be reduced and the strains could be exacted. However, the high-frequency part of the displacement is also eliminated and this is particularly harmful in the case of a high displacement gradient. The measured strains at the boundary also frequently contain relatively large errors compared with those inside the boundary. So it is meaningful to have a method that offers all the advantages of DIC yet is capable of directly measuring stress gradients in the whole field near stress concentrators. Accordingly, as an extension of the DIC-based method, authors proposed a full-field optical technique called Digital Gradient Sensing (DGS) [5,6] which can quantify the elastooptic effect using the DIC method for mechanical characterization of optically transparent planar solids. The method can link angular deflections of light rays quantified using DIC to two orthogonal stress gradients under plane stress conditions. So, the DGS measures Cartesian stress gradient components directly; subsequently, Chandru Periasamy [5] calculated the sum in Cartesian stresses or radial stress in the special case. Significantly, sum stress can be obtained by numerical integration from the stress gradient. Using the integral method has the following advantages compared with the differential method. Differential has an amplification effect on the high frequency signal. For example, it usually needs filtering and smoothing in the process of image processing while integral suppresses the high frequency part, and has no special requirements on the smooth processing. However, individual stress components that directly relate to the material mechanical properties are preferred, and a stress separation method is essential. This is a similar issue as that for photoelasticity, which measures the principal stress difference and the principal direction, and thus the stress components themselves cannot be obtained directly. Several methods have been proposed for stress separation in photoelasticity such as the shear difference method [7,8]; numerical analysis of the compatibility equation, oblique incidence and hybrid methods of photoelasticity, and another experimental method [9–12]. However, all of the above methods for stress separation in photoelasticity are not suitable for that in the DGS method because shear stress is not available when the shear difference method implants in the DGS method. Hybrid experimental methods require a significant modification in the experimental setup, which is a cause of inconvenience, and the other two methods are also impracticable due to their
difficult measurement theory. On the other hand, several hybrid methods of photoelasticity and numerical analysis [13–15], and inverse analysis methods have also been proposed for the stress separation. D. G. Berghaus [13] provides a hybrid method with photoelasticity and finite-element combined, from which displacements along free surfaces and axes of symmetry are obtained from photoelasticity and then perform the FEM to get stress components. S. Yoneyama [15] also provides a hybrid method with photoelasticity and finite-element combined, from which tractions along boundaries are obtained from both principal stress difference and principal direction, and then perform the FEM to get the stress components. In this study, an inverse method is proposed for stress separation in the DGS. First, the stress sum of each point is estimated from the stress sum gradient using the linear least squares method. Then, the interesting area is subdivided into quadrilateral elements, nodal forces along boundary are obtained by the stress sum in the DGS and individual stress components are obtained by the FEM directly. Finally, the inverse method is applied to a three-point bending specimen under compression test. Stress components around the crack tip are obtained; stress components at each point are compared with the FEM direct analysis for this problem; and the results are summarized and conclusions are drawn. 2. Theoretical Background A. Experimental Procedure
Figure 1 schematically illustrates the typical experimental setup for the DGS. The experiment contains a speckle target, a transparent test specimen, and a digital camera. The target plane surface must have a random speckle pattern. The transparent specimen to be tested is placed in front of and parallel to the target plane at a distance, Δ. The target plane is uniformly illuminated using two white light sources. The camera is placed in front of the test specimen at a long distance, L (L ≫ Δ), with its optical axis normal to the specimen surface, imaging the target plane surface through the specimen in the region of interest. The image captured under no-load is
Fig. 1. Experimental setup for DGS. 10 December 2014 / Vol. 53, No. 35 / APPLIED OPTICS
8351
referred to as a reference image or undeformed image. The image captured under loading is referred to as a deformed image. The speckles on the target plane are compared to the undeformed image due to the local change of the thickness and refractive index at a load. So displacements recorded using DIC can be related to angular deflection fields. These angular deflections can further be related to the inplane stress gradient, as discussed in the following section. The measurement principle of the DGS as shown in Fig. 2 is as follows. The thickness of a transparent specimen is B and refractive index is n0 in a no-load state. Cartesian systems, x; y; z and x0 ; y0 ; z0 , are chosen on the specimen and target plane, respectively, so that the z axis is perpendicular to the specimen (or target plane). When the target plane surface is imaged through a transparent specimen, a point, P, on the target plane corresponding to O on the specimen is imaged in the reference state. When imposing a load, the thickness and refractive index of a specimen are all changed due to local stress; thus, a neighbor point, Q, corresponding to O is imaged in the deformed state. That is, light rays, OP, in a reference state corresponding to OQ in a deformed state. The spatial vector, PQ, with known Δ can be related to angular deflections of light rays. On the assumption of isotropic, linearly elastic, and plane stress conditions, the governing equations of the DGS method are expressed as [5] ∂σ xx σ yy δx ; ≈ α CB ∂x Δ δy ∂σ xx σ yy φy ≈ β CB : ∂y Δ φx
(1)
From above, δx and δy can be derived from DIC, Δ is measured in advance, and C and B as material property are also known. So, angular deflections as well as stress gradients are obtained. Stress sums σ xx σ yy , determined from measured stress gradients, are shown in Section 4.
B. Inverse Boundary Method Using Cartesian Stress Components
The interesting region for stress analysis is subdivided into quadrilateral elements as shown in Fig. 3. For a finite element model, if the boundary conditions are applied appropriately, then stress components at every interesting point could be obtained accurately. However, for the local finite element model matching the analysis area, the boundary conditions are usually unknown. The unknown boundary nodal forces around the analysis region are inversely calculated using the measured sums of the Cartesian stress components at the interior interesting points. Then, individual stress components at interior points are determined by finite element direct analysis using the computed boundary conditions. The basic theory of the hybrid method is shown in Fig. 3. The displacements of two orthogonal directions at node, c1 , and the displacements of the y direction at node, c2 , are fixed, so rigid body motion is not allowed. Apply a unit force on an element node along the boundary with its direction in accordance with the x or y axis. Do the finite element procedure and then the stress components at any interesting point can be calculated at this applied unit force. The procedure is repeated by changing the applied node and direction. That is, for the applied load of Pj 1Nj 1…n, stress components of σ 0x ji ; σ 0y ji i 1…m at any interesting data point such as Qj xi; yi are obtained, where j is the unit force index, i is data point index, n is the total number of the unknown nodal forces to be determined, and m is the total number of interesting point. Thus, the stress components of σ x ji and σ y ji at any point xi; yi under real nodal force can be expressed as follows: σ x i σ y i
n X j1 n X j1
σ 0x ji F j i 1…m; σ 0y ji F j i 1…m;
(2)
where F j is the real nodal force along the boundary.
Fig. 2. Schematic of the working principle of the DGS. 8352
APPLIED OPTICS / Vol. 53, No. 35 / 10 December 2014
Fig. 3. Finite element model with unit force applied on the boundary.
So the sum of the stress components (stress sum) obtained by the DGS method can be related to unknown nodal force as follows: σ x σ y i
n X σ 0x ji σ 0y ji F j i 1…m;
(3)
fatigue loading in pure mode I at 1 Hz, average load of 1 kN, and alternating load of 0.3 kN, resulting in a total crack length of 12.2 mm. Then, black paint was sprayed on the target plane surface in order to apply a high contrast feature on the specimen surface.
j1
B. Sum Stress Gradient Measurement
where σ x σ y i is the stress sum obtained by the DGS at the point xi; yi under an unknown applied load of F j. Eq. (3) can be written in matrix form as T bF;
(4)
where 2 6 T6 4 2 6 b6 4
σ x σ y 1 .. .
3 7 7 5
σ x σ y m σ 0x 11 σ 0y 11 .. .
σ 0x 1m σ 0y 1m 3 F1 6 . 7 7 F6 4 .. 5: 2
..
σ 0x n1 σ 0y n1 .. .
.
3 7 7 5
σ 0x nm σ 0y nm
Fn The solution of Eq. (4) in the least squares sense is F bT b−1 bT T:
(5)
Equilibrium equations relate the total applied loads in the horizontal and/or vertical directions to the loads produced by stresses at experimental data points. Using the above equations, the applied loads (nodal forces) along the boundary are calculated. Stress components at every experimental data point are then computed under the calculated nodal forces.
The PMMA specimen is placed in front of a target plane painted with a black and white speckle pattern with (20 mm) distances. A CCD is placed at a large distance (L ≈ 680 mm) from the specimen is combined with a computer for capturing images. At a small load of several Newtons, the image is recorded using a camera (2048 × 2448 pixels) as a reference image. Subsequently, speckle images at different load levels were captured as deformed images. As fields around the crack tip suffer deformation, light rays passing through the specimen were deflected by nonuniform stress, thus the speckle image was deformed relative to the reference state. Displacements were analyzed by the DGS, so full-field displacements and, subsequently, angular deflections (stress gradients) were extracted. During the DGS analysis of these images, an area of 7.69 mm × 5.77 mm around the crack tip is chosen to be the region of interest, as shown in Fig. 4. In the image, 1 mm corresponds to 26 pixels. Angular deflections were calculated with subset size of 29 × 29 pixels and grid step of 5 pixels, which yielded a total number of 1271 with an array of 31 × 41 data points. When calculating stress components we consider rigid body motions and imposing boundary conditions of the problem to quantify the analysis of the contour levels. That is, in the current problem, the boundary conditions, such as asymmetric stress gradients (φy ) in the y direction about the x axis, symmetric stress gradients (φx ) in the x direction about the x axis, and vanishing stress gradients far away from the loading point are all utilized. Thus, Fig. 5 shows the resulting stress gradient of ∂σ xx σ yy ∕∂x and ∂σ xx σ yy ∕∂y with the rigid body translation deleted in the end for the representative load of 600 N in a square region around the crack tip, and the crack tip position is (30.5 mm, 34.5 mm).
3. Experimental Procedure A.
Specimen and Equipment
To verify the proposed hybrid method, a stress concentration problem was studied. An edge-notched three-point bending specimen (220 mm × 60 mm × 9 mm) was manufactured from polymethyl methacrylate (PMMA). The elastic modulus is 2.8 GPa and Poisson’s ratio is 0.37, which are measured by strain gauges; and the stress optical constant is −1.0 × 10−10 m2 ∕N, which calibrated in the Appendix described in the literature [16]. In order to simulate that the crack exists in the actual components, the edge-notched PMMA specimen was subjected to
Fig. 4. Experimental analysis area. 10 December 2014 / Vol. 53, No. 35 / APPLIED OPTICS
8353
Fig. 5. Contours of the sum stress gradient for the x direction (left) and the y direction (right) for 600 N.
C.
2
Sum Stress Calculated by Linear Least Squares
Numerical integration converts gradient data into sum stress by using least squares fitting, pi;j1 pi;j T i;j1 − T i;j ; 2 xi;j1 − xi;j i 1; …; M;
J 1; …; N − 1;
qi1;j qi;j T i1;j − T i;j ; 2 yi1;j − yi;j i 1; …; M − 1;
J 1; …; N;
where x, y, are the Cartesian coordinates, p and q denote the stress gradients in the x and y directions, respectively, i.e., p ∂σ xx σ yy ∕∂x and q ∂σ xx σ yy ∕∂y. T is the stress sum, i.e., T σ xx σ yy . M and N are the total rows and total columns about the calculated area. According to Eq. (6), the stress sum information can be acquired as 1 T i;j1 − T i;j pi;j1 pi;j xi;j1 − yi;j ; 2 i 1; …; M; J 1; …; N − 1; 1 T i1;j − T i;j qi1;j qi;j xi1;j − yi;j ; 2 i 1; …; M − 1; J 1; …; N:
(7)
Furthermore, Eq. (7) can be rewritten in terms of matrices as
6 6 6 T 6 6 6 4 2
where 8354
APPLIED OPTICS / Vol. 53, No. 35 / 10 December 2014
0 −1 1 3 T 1;1 7 T 2;1 7 7 7 ; .. 7 . 7 5
T M;N
;
M−1NMN−1×MN
MN×1
3 p1;2 p1;1 x1;2 −x1;1 6 7 6 p1;3 p1;2 x1;3 −x1;2 7 6 7 6 7 .. 6 7 6 . 7 6 7 6 p 6 M;N pM;N−1 xM;N −xM;N−1 7 7 G 6 7 q2;1 q1;1 y2;1 −y1;1 6 7 6 7 6 7 q y −y q 3;1 2;1 3;1 2;1 6 7 6 7 6 .. 7 6 . 5 4 qM;N qM−1;N yM;N −yM−1;N
:
M−1NMN−1×1
Least squares estimation can be carried out to provide the stress sum as an optimized solution. The counters of the stress sums thus obtained are shown in Fig. 6. D.
AT G;
3
6 0 −1 0 0 1 0 0 7 7 6 6 . . . . . . . . .7 6 . . . . . . . . .7 6 . . . . . . . . .7 7 6 6 0 0 −1 0 0 1 7 7 6 A 6 7 6 −1 1 0 0 7 7 6 6 0 −1 1 0 0 7 7 6 6 . . . . . . . . .7 6 . . . . . . . . .7 4 . . . . . . . . .5 2
(6)
−1 0 0 1 0 0
Stress Components Obtained by Inverse Boundary
The stress separation is performed in the 8 mm × 6 mm region around the crack tip. This area is divided into a finite element mesh which is a little larger than the experimental areas in order to
Fig. 8. Nodal force along a boundary. Fig. 6. Estimated stress sum obtained by the linear least squares method for 600 N.
include all experimental data points. Using the procedures as Section 2 described to calculate the stress for each experimental data point x; y, it must be found to which element it belongs and mapped from global coordinates x; y to local coordinator ε; η. Numerical interpolation is adopted here to obtain the discrete expressions of the mapping [17], considering the fact that it is generally not possible to determine the explicit expressions of the mapping. Figure 7 shows the finite element model of the 8 mm × 6 mm region used for the proposed method. In this model, 4-noded isoparametric elements are used. The numbers of the elements and the nodes are 300 and 336, respectively. In order to obtain the stresses under the unit force at a point on the boundary, the displacements at some nodes must be fixed to prevent the rigid body motion. Displacements of the y direction at node c1 and c2 are fixed; displacements of the x direction at node c3 are fixed. So, rigid body motion is not allowed, and other nodal force on the boundary needs to be calculated. The total number of nodes on the boundary is 70, so the nodal force to be determined is 137 because the three displacement components at the points, c1 , c2 , and c3 are fixed. The nodal forces at 137 points are simultaneously determined as described in Section 2. The nodal forces at 19 points along boundary, c4 ; c5 obtained from the estimated stress sum, are shown in Fig. 8. For comparison, the nodal force obtained from the FEM direct analysis are also shown. As
Fig. 7. Finite element model for analysis area.
shown in this figure, the significant difference between the nodal forces obtained using the proposed method and direct analysis is not observed. Stress components could then be obtained by finite element analysis directly. Obtaining the nodal force and acquisition of stress components are performed using Fortran language programs made by the authors. The numbers of finally separated stress components obtained from the experimental stress sums are shown in Fig. 9. An FEM analysis of σ x and σ y is also shown in Fig. 10. By comparison, we can see that the calculated stress component distribution by experimental stress sum corresponds well to theoretical calculations. E. Discussion of the K -Dominance Characterization About the Selected Area
The optical measurements are usually used in evaluating crack-tip parameters. Crack-tip parameters such as stress intensity factor (SIF) can be determined from full-field in-plane optical data based on the generalized polynomial functions proposed by Williams. A K-dominant field can be defined as the one in which the contribution from the higherorder terms is negligible when compared to the first term. If there is good agreement between experimental and theoretical values when only the first term in the expansion field is used, it suggests that the use of the K-dominant field itself may be sufficient in obtaining an accurate parameter in the specimen’s geometry. However, due to the finite size of the specimens, higher-order nonsingular deformations generally influence the crack-tip fields. Hence, the assumptions of K-dominance are relaxed and higher order terms are included. According to the literature [18,16], there exists a 3D zone surrounding the crack tip and a possible lack of K-dominance. The maximum extent of the 3D region corresponds to approximately half of the specimen thickness. So, the selected area of 7.69 mm × 5.77 mm around the crack tip might lack K-dominance. 3D effects generally prevent data from within r∕B 0.5 (r is the distance from the crack tip and B is the specimen thickness) from being analyzable using 2D 10 December 2014 / Vol. 53, No. 35 / APPLIED OPTICS
8355
Fig. 9. Stress components of σ x and σ y for 600 N determined by the proposed inverse method.
Fig. 10. Stress components of σ x and σ y estimate by finite element direct analysis.
descriptions. For mode I cracks as the mode in this paper, however, both optical and finite-element [19] results have revealed that there exist sectors in the crack-tip vicinity where 3D effects are minimal and the 2D singular term is dominant. The data in these regions, although well within the r∕B 0.5 limit, can be well described by the K-dominant field. So mode I fracture is dependent on the critical values of K I . 4. Conclusions
The inverse method theory for DGS stress analysis is developed using the sum of the Cartesian stresses at interesting data points, which can be calculated from the DGS method when aided by analytical expressions. The DGS measurements of the stress gradient in Cartesian stresses at all interesting points are used to compute the stress sum of each point, and stress sums are used for the inverse boundary analysis. Then the individual stress components are calculated by the FEM directly using the computed boundary conditions. The method has also been implemented to study a stress concentration problem about the crack tip. The results indicate that the actual boundary nodal forces can be predicted with a reasonable accuracy, and stress components can also be obtained reasonably. The work is supported by the National Basic Research Program of China (under Grant 14118683); 8356
APPLIED OPTICS / Vol. 53, No. 35 / 10 December 2014
China Postdoctoral Science Foundation (under Grant 2014M552548XB); and the National Natural Science Foundation of China (under Grant 11402103). References 1. L. B. Meng, G. C. Jin, and X. F. Yao, “Application of iteration and finite element smoothing technique for displacement and strain measurement of digital speckle correlation,” Opt. Lasers Eng. 45, 57–63 (2007). 2. M. A. Sutton, J. L. Turner, H. A. Bruck, and T. A. Chae, “Fullfield representation of discretely sampled surface deformation for displacement and strain analysis,” Exp. Mech. 31, 168–177 (1991). 3. D. J. Segalman, D. B. Woyak, and R. E. Rowlands, “Smooth splike-like finite-element differentiation of fullfield experimental data over arbitrary geometry,” Exp. Mech. 19, 429–437 (1979). 4. B. Pan, H. M. Xie, Z. Q. Guo, and T. Hua, “Full-field strain measurement using a two- dimensional Savitzky–Golay digital differentiator in digital image correlation,” Opt. Eng. 46, 033601 (2007). 5. C. Periasamy and H. V. Tippur, “A full-field digital gradient sensing method for evaluating stress gradients in transparent solids,” Appl. Opt. 51, 2088–2097 (2012). 6. C. Periasamy and H. V. Tippur, “Measurement of crack-tip and punch-tip transient deformations and stress intensity factors using Digital Gradient Sensing technique,” Eng. Fract. Mech. 98, 185–199 (2013). 7. S. J. Haake and E. A. Patterson, “The determination of principal stresses from photoelastic data,” Strain 28, 153–158 (1992). 8. M. Solaguren-Beascoa Fernández, J. M. Alegre Calderon, P. M. Bravo, and I. I. Cuesta Segura, “Stress-separation techniques in photoelasticity: a review,” J. Strain Anal. Eng. Des. 45, 1–17, (2010).
9. R. J. Greene and E. A. Patterson, “An integrated approach to the separation of principal surface stresses using combined thermophoto-elasticity,” Exp. Mech. 46, 19–29 (2006). 10. T. Sakagami, S. Kubo, and Y. Fujinami, “Full-field stress separation using thermoelasticity and photoelasticity and its application to fracture mechanics,” JSME Int. J. Ser. A 47(3), 298–304 (2004). 11. S. Yoneyama, Y. Morimoto, and M. Kawamura, “Twodimensional stress separation using phase-stepping interferometric photoelasticity,” Meas. Sci. Technol. 16, 1329–1334 (2005). 12. J. Lim and K. Ravi-Chandar, “Dynamic measurement of two dimensional stress components in birefringent materials,” Exp. Mech. 49, 403–416 (2009). 13. D. G. Berghaus, “Combining photoelasticity and finite element methods for stress analysis using least squares,” Exp. Mech. 31, 36–41 (1991). 14. D. Chen, A. A. Becker, I. A. Jones, T. H. Hyde, and P. Wang, “Development of new inverse boundary element techniques
15. 16. 17. 18.
19.
in photoelasticity,” J. Strain Anal. Eng. Des. 36, 253–264 (2001). S. Yoneyama, S. Arikawa, and Y. Kobayashi, “Linear and nonlinear algorithms for stress separation in photoelasticity,” Exp. Mech. 52, 529–538 (2012). R. Zhang, R. Guo, and S. Wang, “Mixed mode fracture study of PMMA using digital gradient sensing method,” Eng. Fract. Mech. 119, 164–172 (2014). Y. P. Yu, Y. Xiang, Q. S. Liu, and A. Zhong, “Research on numerical inverse isoparametric mapping interpolation and its application,” Rock Soil Mech. 22, 226–228 (2001) (in Chinese). S. Ramaswamy, H. V. Tippur, and L. Xu, “Mixed-mode cracktip deformation studied using a modified flexural specimen and coherent gradient sensing,” Exp. Mech. 33, 218–227 (1993). H. V. Tippur, S. Krishnaswamy, and A. J. Rosakis, “Optical mapping of crack tip deformations using the methods of transmission and reflection coherent gradient sensing: a study of crack tip K-dominance,” Int. J. Fract. 52, 91–117 (1991).
10 December 2014 / Vol. 53, No. 35 / APPLIED OPTICS
8357
