Inverse problems with the digital image correlation: approach and applications

problems in It on regression algorithms to estimate materials and/or loading parameters experimentally-evaluated displacement


INTRODUCTION
ptical techniques have been used over a time-scale of more than one century as powerful non-contact tools to determine the mechanical properties of materials. Photoelasticity [1], for instance, was employed for a variety of stress analyses and even in routine design approaches especially before the advent of numerical methods. It has been proven to be particularly useful to investigate highly localized stress state [2][3][4][5][6][7][8]. Interferometric approaches were exploited to measure the displacement field experienced by samples and over-deterministic methodologies, aimed at the calculation of the elastic constant of brittle materials, were also proposed [9,10]. Unfortunately, the application of such techniques, for in-plane stress and displacement/strain measurements, needs special equipment, longtime sample preparation and stringent stability requirements [1]. In addition, they are limited to the use of light-polarizing materials. In the last years, a new full field measurement technique, known as digital image correlation (DIC), was introduced by Sutton et al. [11][12][13][14][15][16][17][18][19] and it has been widely used in supporting the characterization of engineering alloys. This technique offers attractive advantages as it can be applied to any kind of material; the only requirement is a surface pattern for pixels tracking which can be easily obtained by airbrushing, or simply exploiting the natural texture of the samples. Thanks to these interesting features, DIC has been used for different applications to analyze the strain field in components during operating conditions [20][21][22], to characterize the strain/stress level in biomedical implants [23,24] or to estimate the elastic properties of materials [25][26][27]. Moreover, DIC was used to calculate the fracture properties of materials [28][29][30] starting from the crack tip displacements. A procedure to calculate the mode I fracture parameters for an orthotropic body starting from full-field measurements provided by DIC was applied in [31,32]. The mixed mode generalized stress intensity factors (GSIF) was experimentally evaluated from DIC displacement and strain fields by means of a path independent integral in [33]. An over-deterministic approach, based on the least square regression of the DIC measured displacement, was also applied to estimate the stress intensity factor of non-conventional materials such as Shape Memory Alloys (SMA) and super-alloys [34][35][36][37][38][39][40][41][42][43], where common elastic and/or elastic-plastic theories cannot be directly applied. Similar approaches have been used for the estimation of the elastic constants of brittle materials by means of Brazilian tests [44][45][46][47]. However, it has to be pointed out that most of these works are based on a simple linear regression that fails when plasticityinduced non linearity mechanisms are involved. In this paper a reliable approach is proposed to solve inverse problems using optical techniques, even accounting for material non-linearity. The methodology is based on DIC measurement of the displacement field experienced by a sample during a thermo-mechanical test and the application of the regression methods to calculate unknown parameters such as material constants and/or external loads. An automatic procedure, based on displacement data, was developed and validated by means of finite element simulations. It is able to estimate the parameters of interest accounting for the unavoidable rigid body motions and to fit the experimental displacement fields to representative analytical solutions. In addition, thanks to iterative calculation algorithms plasticity-induced non-linear effects, can be captured successfully minimizing the estimation errors. This latter represents an important aspect of the methodology as it allows to account for important mechanisms in structural mechanics that are impossible to capture with more traditional DIC-based methods. An example is represented by the plasticity-induced crack opening/closure phenomena in fracture mechanics problems. The knowledge of the effective crack tip position, in fact, is fundamental [48] for a correct estimation of the fracture properties of materials. Moreover, the reference system position of the displacement field, referred to the analytical solutions, is another crucial parameter which need to be properly identified within the recorded images for the application of the regression approach. Very often, it is complex to be visually captured because of geometry and/or resolution-related constraints. This is particularly evident for axisymmetric problems. The ad-hoc developed iterative procedure allows to estimate the reference system position accurately and automatically. The methodology was applied to three case studies with the aim to evaluate: i) the stress intensity factor on fracture mechanics problems, ii) the elastic properties of a material by means of Brazilian test, iii) the contact pressure generated by thermally activated shape memory alloy (SMA) rings used for pipe coupling [49]. Direct comparison between the regressed solutions with those based on conventional techniques in experimental mechanics confirmed the accuracy of the proposed methodology. Results revealed that this latter represents a viable technique for materials/components characterization.

ESTIMATION OF THE UNKNOWN PARAMETERS: LEAST SQUARE METHOD
he determination of the parameters of interests can be performed by the regression analysis of the displacement fields measured by the DIC. For a 2D elastic problem, the displacement experienced by a body in a certain point can be expressed as follows: It is worth noting that the vector {U} can represent either the material properties and/or the loading conditions. Therefore, depending on the specific requirements, the method can be used for material characterization or for estimating the applied thermo-mechanical loads.

DISPLACEMENT FIELDS
n next paragraphs, the analytical relations of the displacements, for all the case studies, are reported and described with respect to the most general formulation reported in Eqn. (1).

Case study 1: near crack tip displacement field
Williams' solution, for an elastic homogeneous 2D cracked body [50], provides an approximation of stress and deformation fields by means of their expansion into power series. According to the Eqn. (1) with n=∞, it follows that: If the only first term of the Williams' series is considered (n=1), in the case of mode I loading, Eqn. (4) becomes: where KI represents the mode I stress intensity factor (SIF) and [] terms are the following: It is worth noting that these equations are suitable in describing the displacement field in a very close crack tip region, i.e. the KI dominant zone [48]. If a larger area is investigated, the analytical description of the displacement field requires the introduction of further terms. Thanks to the high number of data recorded by the correlation technique, the methodology can provide not only the SIF, which represents the most common fracture parameters, but also the non-singular higher order terms called T-stress. Most recently, it was shown that this latter can also play a significant role in the analysis of fracture behavior. In fact, the sign and magnitude of the T-stress modifies the size and shape of the crack tip plastic zone [51,52] and affect the stability of the crack path direction [53][54][55]. Hence, the application of both the SIFs and the T-stress in failure investigations is becoming increasingly popular in the scientific community and many efforts are being made to calculate their values. Therefore, if a further term, n=2, is added to the Eqn. (8), the new equations can be written as follows: The theoretical displacement fields, calculated from Eqn. (10), are illustrated in Fig. 1 in the form of contour maps.

Case study 2: displacement fields of a brazilian disk
The displacements field of a disk, under plane stress condition, subjected to a compression load P along the y direction (see Fig. 2), can be expressed as follows [56]: In Eqn. (12) coordinates   , x y are referred to the center of the disk (x 0 , y 0 ), t is the disk thickness, D is the disk diameter, (r 1 ,  1 ) and (r 2 ,  2 ) are the polar coordinates of the generic point with respect to the contact points A and B reported in Fig.   2. These latter can be expressed as follows: The theoretical displacement fields, calculated from Eqn. (12), are illustrated in Fig. 2 in the form of contour maps. Figure 2: Schematic depiction of a Brazilian disk subjected to a compression load P, together with the corresponding displacements along the x and y axis.

Case study 3: displacement field of an axisymmetric component subjected to pressure and thermal load
The displacement field of an axisymmetric component subjected to an external pressure P e , see Fig. 3, can be written as follows [56]:  The total displacement utot, can be obtained by summing the two displacement components, utot=√[(ux) 2 +(uy) 2 ]. Its trend is reported in Fig. 3. If the system is also subjected to a thermal load, a new term has to be added to the Eqn. (18).
cos sin r r (21)  is the thermal expansion coefficient and T is the temperature variation. It is important to observe that, similarly to the previous case, the methodology can be used to accomplish two purposes: i) estimation of the applied load P e and T if the material under investigation is known, ii) estimation of the properties of the material, E,  and , if the applied loads are measured by proper transducers. In this work the unknown parameters are P e and .

SENSITIVITY ANALYSIS
qn. (2) and (3) shows that the points of interest must be selected in such a way that the obtained equations are linearly independent, avoiding ill-conditioning of the pseudo inverse of the matrix To limit such mathematical issue, an over-deterministic approach exploiting a large number of data points can be used. In fact, considering the displacements of the whole domain, many equations as much are the investigated points can be considered and estimation errors of the unknown parameters are significantly reduced. However, it is important to identify the displacement components to be used for minimizing the estimation errors. To this aim, a sensitivity analysis, was performed using the following cost functions, 2 Φ x u and 2 Φ y u : , , In conclusion, uy displacements are recommended for the fracture parameters calculation (case study 1) while ux displacements are preferred for the estimation of the elastic constants (case study 2). Both displacements can be used for contact pressure and thermal expansion coefficient calculation (case study 3). However, for more accurate results, both the fields can be considered.

OF THE REFERENCE SYSTEM
qns. (10,12,20) represent the solution of the displacements when the coordinates (x 0 , y 0 ) of the reference system are known and no rigid body motion are involved. When dealing with experiments, such conditions cannot be satisfied because the experimental displacements, provided by DIC commercial software, are typically referred to a system different from the required one (x 0 , y 0 ) and because unavoidable rigid body motions always occur. Therefore, they have to be included in the calculation process. Rigid body motions can be easily added in the most general equation of the displacements, Eqn. (1), as follows: 11  12  13  14  15  1  2  21  22  23  24 25 where A is the rigid body rotation, B x and B y represent the rigid translation along the x and y axis, respectively. The addition of the latter parameters does not represent an issue for implementing the over-deterministic method because the displacements {u} can be still expressed by a linear combination of the unknown parameters {U} and the number of data, provided by the correlation technique, is still significantly higher than the number of unknowns. The estimation of the reference frame location (x0, y0), instead, cannot be simply performed because the introduction of such coordinates within the vector of unknowns {U}={U 1 U 2 A B x B y x 0 y 0 } T leads to a system of non-linear equations to be solved by a non-linear fitting procedure based on Newton-Raphson method. To this aim, eqs. (10,12,20) can be written as a series of iterative equations based on Taylor's series expansions as follows: where the matrix [] i (2 × 7) represents the gradient of the displacement vector {u} i (see Eqn. 26). If Eqn. (27) is applied to the m measurements points and the displacement vector at the i+1 step is set to the experimental one ({u * } i ={u * } -{u} i ) the following overestimated system of 2m equations is obtained: The procedure described above is repeated until the corrections {U} become acceptably small. It is important to point out that convergence cannot be easily reached if the initial trial values for {U} are too different from the actual ones.

NUMERICAL VALIDATION
n order to validate the proposed approach, finite element simulations were carried out. 2D linear elastic problems were analyzed. Four-nodes shell elements were used and special care was done in order to guarantee a regular mesh on the region of interest (selected domain used for the proposed numerical procedure). Fig. 7 shows a depiction of the models together with the main dimensions and boundary conditions. All the geometries are one millimeter thick. The material properties are listed in Table 1. The least squares regression was performed on the FEM displacements and Eqn. (24) was used to calculate the unknown parameters Ui. For the Brazilian disk and the ring, results were compared with data, reported in Table 1, used as input data for the simulation, i.e. E=5 GPa and =0.4 for the disk and P e =150 MPa and =1.3‧10 -5 °C -1 for the ring. For the single I edge crack sample, the regressed value of the SIF was compared with that calculated using the standard formula (see Eqn. 30) based on a linear elastic fracture mechanics approach [48].  Table 1: Material properties and loading conditions used for the numerical simulations. Figure 8: Comparison between the numerical and regressed displacements: a) ux displacements for cases study 1; b) uy displacements for cases study 1; c) u x displacements for cases study 2; d) u y displacements for cases study 2; e) u tot displacements for cases study 3.
Figs. 8a-8e show a comparison between the numerical displacement contour plots (blue lines) and the regressed ones (red line). Results show good agreement and the regressed parameters, Ui, revealed an error lower than 1% compared to the reference values. Moreover, it is important to note that the methodology was able to automatically identify the location of the reference system for each case study.

MATERIALS AND EXPERIMENTS
he displacement field experienced by the samples, in all the case studies, was in-situ monitored during thermomechanical tests by using a CCD Camera (Sony ICX 625 -Prosilica GT 2450) with a resolution of 2448x2050 pixels. In addition, a suitable objective was adopted to focus the region of interest (Rodagon f. 80 mm -Rodenstock). Digital Image Correlation was performed by using a commercial software (VIC-2D®, Correlated Solutions). T Next sections describe the material samples and testing conditions used in the investigated case studies.

Case study 1
The first case study consists in evaluating the SIF and/or fracture toughness of a sample. To this aim, Single Edge Crack (SEC) specimens, with dimensions shown in Fig. 9, were manufactured from aluminum AA 6005-T6 sheets with thickness t =1.5 mm, by Electro Discharge Machining (EDM). The rolling direction was parallel to the tensile axis. The samples were fatigue pre-cracked (f=5 Hz, R=min/max=0 and max=20 MPa), starting from EDM notch (radius around 100 μm), up to a length to width ratio, a/W, close to 0.40. Almost straight crack path, normal to the load direction initiating from EDM notches, was obtained. Fig. 9 also reports an image of the real sample together with the intensity histogram of the speckle pattern. Tensile loading-unloading cycles were performed by using an electro-dynamic testing machine (Instron E10000), applying a maximum load P=900 N, corresponding to a maximum stress max=Pmax/Wt=50 MPa.
Correlation analyses were performed by using a subset size of 41 pixels and a subset distance of 5 pixels. Finally, isothermal displacement-controlled (0.05 mm/min) fracture tests were also carried out by monotonic tensile loading until fracture. In order to prove the reliability of the measurements, a linear elastic fracture mechanics approach [48] was adopted to calculate the SIF by using the following equation:  In addition, monotonic tensile tests were performed on dog bone shaped samples to characterize the elastic properties of the material (E=73 GPa, y=295 MPa), necessary for the application of the methodology. They were manufactured from the same aluminum sheets. Three samples, with 1.5 mm of thickness and 3.5 mm of width, according to ASTM E8, were used to perform the tests. Displacement-controlled (1 mm/min) mechanical tests were carried out by an electro-dynamic testing machine (Instron E10000) equipped with a 10 kN load cell and an extensometer with gauge length of 10 mm was used to measure the deformations.

Case study 2
In the second case study, the Brazilian disk test was combined with the proposed DIC-based procedure to calculate the elastic constants, i.e. the elastic modulus and the Poisson's ratio, of a PVC polymer. Brazilian test is a valid test to evaluate the mechanical properties of brittle materials due to the well-known difficulties to perform more common tensile tests.
An electro-mechanical testing machine (MTS Criterion s42, USA) equipped with a 5 kN load cell has been used. All the experiments were carried out in displacement control mode (0.05 mm/min) and curved loading platens with 42.5 mm of radius, according to the ISRM [57], were used (Fig. 10). Six polymeric samples were investigated for statistical aim and all of them were airbrushed with white and black paint in order to get a proper pattern with a random grey scale distribution. Fig. 10 reports an image of the speckle pattern of the sample together with the intensity histogram. Samples were made with a diameter of 55 mm ad a thickness of 8 mm in order to get plane stress conditions, according to the Timoshenko theory [56]. Correlation analyses were performed by using a subset size of 61 pixels and a subset distance of 5 pixels.
In addition, tensile tests were performed on the same material in order to calculate the effective elastic properties; three dog-bone shaped samples with 8 mm of thickness and 19 mm of width, according to ASTM D638, were used to perform the tests. HBM K-XY3 strain gauges were glued at the center of the samples and the HBM QuantumX MX1615B Strain Gauge Bridge Amplifier together with the Catman software was employed to record longitudinal and transverse strains. The crosshead speed was set to 1 mm/min and the Poisson ratio, ν, was calculated taking into account the transverse sensitivity of the strain gauges and the Poisson ratio of the material adopted for their calibration.

Case study 3
The last case study consists in the measurement of the coupling pressure generated by a SMA ring when it is thermally activated on axisymmetric components. SMA rings exploit temperature-induced phase transformation mechanisms to offer new ways to join and seal cylindrical components as demonstrated by some recent research and/or technological studies as well as by a significant number of patents [58][59][60][61][62]. A 38NiCrMo3 steel ring, simulating the pipe system, was clamped with an external commercial Ni46.5Ti45.0Nb8.5 SMA ring. These metal rings are supplied in the expanded condition. Once shrunk, they apply a uniform gripping pressure that is onto the substrate. Thanks to the thermal hysteresis of the material, the radial force is kept during the subsequent cooling of the assembly to room temperature, which represent the operating condition of the coupler. Fig. 11 illustrates the geometry and dimension of the NiTiNb SMA ring used in this investigation, together with a schematic depiction of the SMA-steel ring assembly analyzed for contact pressure measurement. Figure also reports an image of the speckle pattern of the steel ring together with the intensity histogram. An initial SMA ring/steel rings diametric clearance of about 0.15 mm has been chosen. In order to prove the accuracy of the least square method in calculating the radial force generated by the SMA ring on the steel one, strain gauge (SG) measurements were also performed. In particular, x-y electrical strain gauges (XC11, HBM), having a wide operative temperature range (-200 °C/+250 °C), were mounted at the internal diameter of the steel ring (see Fig. 11) by using a two component epoxy resin (EP310S, HBM). Two strain gauges were glued at a relative angle of 120 °C (see Fig. 11) and quarterbridge configurations were adopted, which gives four different signals from the two x-y SGs. In order to delete the apparent thermal strain related to the applied temperature variation, namely the thermal output, preliminary experiments using the only steel ring, i.e. with no applied contact pressure, were performed. The wall thickness of the steel ring was chosen with the aim of maximizing both the strain measurements sensitivity, obtained by electrical strain gauges (see Fig. 11), and the displacements sensitivity for the digital image correlation measurements.
The coupling system was heated up to T=200 °C and then cooled down at room temperature (T op ) in order to guarantee a proper SMA ring activation [59]. A k-type thermocouple was used to monitor the steel ring temperature. Correlation analyses were performed by using a subset size of 69 pixels and a subset distance of 5 pixels.
SG measurements provided the tangential and axial strain components (and z) at the inner radius of the steel ring and, based on theory of elasticity, the contact pressure P e has been calculated as:

RESULTS
he accuracy of the proposed methodology was demonstrated by systematic comparison of the investigated parameters (material constants or loading parameters) with those obtained by using conventional experimental mechanics techniques.
Case study 1 Fig. 12 shows the evolution of the SIF, K I , as a function of the applied load during a complete loading-unloading cycle. In particular, figure reports both  The measured displacements near the crack tip (blue contours) and those obtained by the regression approach are compared in Fig. 13. In particular, Fig. 13a reports the displacement contours along the x axis whereas Fig. 13b reports the ones along the y axis recorded at P=900 N for a crack length a=3.5 mm. Results show very good agreement. a ) b ) Figure 13: Comparison between the DIC recorded displacement contours (blue curves) and the regressed data (red dashed curves) for a SEC specimen (a=3.5 mm) at P=900 N: a) displacement u x and b) displacement u y . Fig. 14a illustrates load-displacement (P-) curve obtained from the fracture test performed on the aluminum SEC sample. The values of the mode I SIF (K I ), according to Eqn. (30), can be also obtained from the right axis of the figure, as the ratio K I /P is a constant depending on the a/W ratio of the specimen. Figure shows that P- curve exhibits almost linear trend in the early stage of the loading, but evident nonlinear deviations occur near the maximum load (P max ). Due to the high non linearity, involving aluminum alloy when high load is applied, calculation was stopped at the load PQ, measured according to the standard ASTM E399 (see Fig. 14a). In fact, below such load value, the material non linearity is always within the acceptability range of the standard ASTM E399 as shown in Fig. 14a. However, it is important to underline that the aim of this paper is not to evaluate the toughness of the aluminum, but to show the viability and potentiality of the proposed approach even when non-linear mechanisms are involved.
In fact, as shown in Figs. 14, for load values lower than about 1400 N the P- curve is perfectly linear and the regressed and the analytical SIF values are very similar. For higher values of the load the P- curve starts to deviate from the linearity as well as the regressed K I . In fact, within this load range, regressed SIFs show higher values than the analytical one. This behavior is a consequence of the non-linear mechanisms occurring at the crack tip when high loads are applied. In fact, due to the high stress concentration, local plasticity occurs generating a redistribution of the stress and displacement field, according to the Irwin criteria [48]. In particular, in standard elastic-plastic metals, the vertical asymptote of the crack-tip stress is shifted by a quantity a, that is proportional to the extent of the plastic zone (a α (KI/Sy) 2 ); this results in an effective crack-tip length and SIF, namely a e and K Ie , as described by the Irwin correction of the LEFM. In this case, the introduction of the parameters x 0 and y 0 in the calculations (cfr. Eqn. 13) is crucial because it allows to find the effective crack tip location and to properly fit the linear region of the experimental displacements, even in case of local plastic phenomena, see Fig. 15.   Fig. 15b, instead, shows the same data when the regression is performed according to eqs. (26)(27)(28)(29). In the latter case, the effective crack tip location was properly identified and the plasticity-induced estimation errors are largely reduced.

Case study 2
Figs. 16 report the comparison between the experimental (blue contours) and regressed displacements (red contours) after rigid body motions contributions (A, B x and B y ) were subtracted. In particular, Fig. 16.a report the u x displacements, Fig.  16.b reports the uy ones. The elastic constants obtained by least square regression, together with those obtained from the tensile tests, are reported in table 2.   A limited mismatch was observed between the results obtained using the regression method and the tensile tests. In addition, a good estimation of the effective position of the disk center was obtained.

Case study 3
Figs. 17a reports a comparison between the total displacements, u tot , experimentally measured by the DIC (blue contours) and the regressed solution (red contours) after the elimination of the rigid body motions, A, Bx and By, and the thermal contribution ‧T‧r of the internal steel ring (=1.3‧10 -5 K -1 ).
In addition, Fig. 17b reports a comparison between the Pe -T curves obtained from strain gauge measurements (Eqn. 31) and from the regression approach (Eqn. 20). A satisfactory agreement, in terms of both displacement and pressure profile, between the two methods was observed.
A maximum difference (around 10%) between DIC and SG measurements in terms of contact pressure estimation has been recorded in the heating stage but with very small differences (never greater than 2%) at room temperature, after the cooling. The higher mismatch observed at high temperatures are probably due to the temperature gradients within the material, i.e. between the bulk and the ring surface where DIC is applied.

CONCLUSIONS
viable approach based on inverse methods was proposed to characterize elastic materials/components or to estimate external loads. The methodology is based on the linear/non-linear regression of the displacement field experienced by a sample subjected to thermo-mechanical loadings. The digital image correlation technique was used to measure the displacement data which were analyzed by ad-hoc developed algorithms to calculate the unknown parameters. Three different case studies were investigated to estimate: i) the stress intensity factor in fracture mechanics problems, ii) the elastic properties of materials by means of Brazilian test, iii) the contact pressure generated by thermally activated shape memory alloy (SMA) rings used for pipe coupling. The proposed procedure was able to estimate successfully the parameters of interest and to fit properly the displacement fields experimentally evaluated. In addition, thanks to ad-hoc developed iterative numerical algorithms, non-linear phenomena such as plasticity can be taken into account in the analyses, minimizing the estimation errors. The unknown parameters calculated by the proposed method were compared with those measured by using conventional experimental techniques and very good agreement was observed, confirming the reliability of the approach.