Computational measurements of stress fields from digital images

Since the pioneering works on digital image correlation, a huge effort has been devoted to the elaboration of strategies coupling measurements with numerical simulations. The main issue achieved with such coupling is the identification of constitutive law parameters from experimental configurations yielding heterogeneous strain/stress fields. The main limitation of this framework is that a constitutive law is to be postulated, what sometimes reveals far from being straightforward. To circumvent this limitation, a methodology is proposed to build a parametric description not only of the displacement field for digital image correlation but also of the local stress. With no assumption on the constitutive relation, the method allows for estimating displacement, strain, and stress fields. The ability of the method is illustrated for different situations.


INTRODUCTION
Since the pioneering work of Lucas, Kanade, and others, 1 image registration or digital image correlation (DIC) has been widely used and developed in many ways. For mechanical experimentalists, the purpose of DIC is to measure displacement fields by comparing digital images of the sample under 2 different loading states. Subset-based techniques are the most popular ones (see, for example, Sutton et al 2 ) because they are the techniques implemented in commercial softwares. One should also mention the non-parametric formulation proposed by Vercauteren et al 3 that allows for measuring a displacement field pixelwise. Compared to the conventional experimental methodology for constitutive behaviour analysis, which consists in performing uniaxial homogeneous stress/strain tests, DIC has the major advantage to provide for full displacement and strain fields. While this is useful to assess the quality of conventional tests, making good use of the rich data extracted from DIC measurements is still the purpose of many research works. One of the main issues is that the stress field in heterogeneous state configurations cannot be accessed without coupling measurements with numerical simulations.
Recently, Besnard et al proposed a finite element (FE) formulation for DIC. 4 This considerably increases the potential of DIC as a strong coupling with FE simulations is established. For example, it is shown in Réthoré et al 5 that both DIC and elastic balance of momentum can be solved at the same time. For the characterization of material constitutive relations, this opportunity to couple DIC measurements with FE simulations has been exploited to design robust identification techniques. 6,7 By identification, we mean here estimating the parameters of a presupposed constitutive model. In some cases, making explicitly a constitutive model and formalizing a constitutive law turn out to be extremely difficult. This is a case, for example, for ductile tearing, as the stress/strain and triaxiality levels occurring at crack tip cannot be reached in homogeneous stress/strain test configurations. At this point, the advantage of this framework (FE DIC combined with FE simulation) reveals also to be a weakness, as it requires the constitutive relation to be explicit.
To circumvent the difficulty of formulating a constitutive equations, manifold learning and dimensionality reduction 8 are good candidates. The idea is to use databases of admissible material state (at least stress and strain) points to extract manifold and then to perform multidimensional interpolation/regression for calculating internal forces in FE simulations. The data-driven computational mechanics framework recently proposed in Kirchdoerfer and Ortiz 9 goes one step further. Instead of minimizing the energy under the constraint that the constitutive relation is satisfied, the authors propose in this work to minimize the distance between the computed material states and those given by a database under the constraint of equilibrium. These non-parametric approaches of the mechanics of materials have a great potential to circumvent the limitations of parametric approaches concerning the formulation of constitutive laws. However, if the database uses results from conventional experiments using homogeneous stress/strain configurations, these promising methodologies would suffer from the same limitations concerning the lack of intense and/or complex loadings.
Full-field measurements could be used on heterogeneous state configuration to enrich material state databases with more complex configurations. However, as it has been mentioned above, the stress analysis of complex tests requires the coupling with FE simulations that need a constitutive model. This is disappointing to observe that these 2 advanced methods (one that extracts material data from experiments and the other that makes use of material state data to perform simulations) cannot help each other. At this point, one ends up with the following fundamental question: Is it possible to measure stress fields? In general, it has not yet been possible. There exist some optical techniques, like photoelasticity, caustics, or more recently diffraction contrast tomography, that allow to measure elastic stress fields, but non-elastic stress evaluation is out of their scope. The aim of this paper is a proposition that may answer this question. It is not claimed that stress may become observable, but we propose herein to formulate auxiliary problems involving DIC so that strain and stress fields can be directly estimated with no assumption on the constitutive relation. Plus, a supplementary information is obtained that can be interpreted as an internal variable field that measures the distance between the actual local constitutive equation and linear elasticity according to a given but arbitrary elastic tensor. The method will be referred to as mechanical image correlation (MIC), as it allows to measure from digital images the local mechanical state of a material.
In the next section, the principle of DIC is briefly explained. Then, the details of MIC are provided in Section 3. In the last section, some examples assessing the feasibility of the method and illustrating its potential are presented.

DIGITAL IMAGE CORRELATION
The DIC framework used in the following is detailed in this section in a similar way as in Réthoré. 7 This formalism is generic but convenient for the development presented in the next section.

Formulation
From 2 grey-level images f and g of a given sample, DIC consists in a non-linear minimization process that identifies the displacement field u that produces the advection of the local texture: where p is the index of a pixel within the region of interest (ROI) over which the computation is performed and x p the coordinates of pixel p on the image frame. For solving this minimization problem, 2 main routes can be distinguished. A non-parametric route is proposed in Vercauteren et al. 3 In this approach, no assumption is made on the variation of displacement and the problem is solved pixelwise using a staggered algorithm with relaxation and smoothing steps. The second route consists in a parametrization of the displacement as follows: In this decomposition,  is the set of degrees of freedom (DOFs) u k , and N k are the shape functions. Different parametrizations can be used. The particular case of crack tip fields required the use of Williams' series as basis functions. 10 Other specific bases, like for 4-point bending test 11 or Brazilian test, 12 can be used, but for arbitrary displacement, FE-based 4,13 or CAD (Computer-Aided Design)-based kinematics 14,15 are more appropriate.
From a first guess for the displacement, an iterative procedure is run until convergence is reached. Assuming that the solution increment du is small enough, a first-order Taylor expansion of g(x p + u(x p ) + du(x p )) is used to linearize the functional. After some manipulations, we obtain that the solution increment is calculated through with and ∇ denoting spatial derivation.
[N] is a matrix that collects the value of the basis functions at the pixels, [∇G] a diagonal matrix that collects the values of ∇g(x p + u(x p )), and {F − G} a vector containing the values of f(x p ) − g(x p + u(x p )). For the sake of clarity, indices indicating the iteration number have been omitted. For a detailed description of the algorithm in a multiscale framework and its implementation, the interested reader may refer to Réthoré et al. 16 To reduce computational cost, a modified algorithm is adopted: We substitute ∇G, which depends on the current solution vector, by ∇F.
The modified system to solve thus becomes with and The advantage of such a modification is that M and parts of b can be computed once for all. Note that when convergence is reached, ∇F equals ∇G. After a few iterations (2 or 3), the modified algorithm and the initial one reveal close to identical.

Condensation
At this point, it is interesting to draw the attention on the following assessment: For any parametrization of the displacement field through N k shape functions, DIC allows for estimating directly from the images the corresponding DOFs u k . If Williams' series are used as in Roux and Hild, 10 a direct access to stress intensity factors is provided. An other interesting example is described in Leclerc et al 6 and Réthoré 7 : The FE displacement is considered as depending on the parameters of the constitutive law of the material. The DIC problem is then solved using not the full FE parametrization but a combination of the FE shape functions corresponding to the sensitivity of the displacement fields to the material parameters. In this context, DIC provides for a measurement of the material parameters. In Section 3, a parametrization of the displacement fields in terms of local stress is described. In other words, using appropriate combinations of FE displacement DOFs, one accesses to the local stress field, which can thus be estimated by solving DIC with the proposed basis of nodal displacement vectors. The general solution algorithm for accessing the estimation of these generalized DOFs by DIC is derived below.
A general condensation method is proposed from an FE parametrization of the displacement. In this case, N k are FE shape functions and u k nodal displacements. Let us consider that the meaningful quantities are not the nodal displacements but linear combination of them. Given a linear condensation operator L, the nodal displacement U is then derived according to In this formalism, collects the value of the generalized DOFs corresponding to the amplitude of the linear combination of nodal displacements defined by L. The DIC problem can be solved using this condensation with the following system: with

MECHANICAL IMAGE CORRELATION
The approaches in Leclerc et al 6 and Réthoré 7 are parametric identification techniques: They allow for estimating experimentally the parameters of a constitutive law from digital images and load measurements corresponding to heterogeneous stress/strain states. These methods give access to a stress field compatible with the presumed constitutive law. As evidenced in Section 1, this is limiting in those cases when explicitly writing the constitutive equations is not possible. The purpose of this section is to propose a formulation giving access to local stress through the parametrization of an anelastic strain field.

Formulation
Let The aim is now to find another route to define these gaps. The actual solution could then be retrieved from the reference solution and the gaps. Because both σ o and σ balance the external load,̄is a self-balanced stress field. It admits the following properties: div̄= 0 over the ROI,̄n = 0 over the boundary of ROI, However, there is a priori no straightforward relation between̄andū. Without loss of generality, one introduces an anelastic strain field a such that̄= wherēis the symmetric gradient of the kinematic gapū. Then using Equation 14, one obtains div( ō) − div( o a ) = 0 over the ROI, ( o (̄− a ))n = 0 over the boundary of ROI.
It is thus clear that if a is known, then the gap from the reference solution to the actual one can be retrieved and the latter can thus be computed. We will then proceed in 2 steps: First, formulate auxiliary problems governing the anelastic strain a , and second, find a parametrization of a so that the auxiliary problems can be solved and the actual solution recovered. Before these 2 steps are described, it is important to point that it is possible to separate a fields in 2 contributions: Once the problem defined by Equation 16 is solved, a solutionū is obtained up to a rigid body motion. It is then possible to define the compatible part I a of a as the symmetric gradient ofū: The remaining part II a = a − I a is not compatible in general, and it is not associated to any displacement (ū The resulting self-balanced stress field is then written as Conversely, I a induces a displacementū I a but no stress as I a =̄. One ends up with 2 families of anelastic strain fields that are able to close the kinematic and static gap, respectively. The following comments arise: • I a being compatible, it is intrinsically related to a displacement (namely,ū). If a set of basis functions deriving from a set displacement functions is chosen for I a , they can be used as a basis for u by using Equation 12.

Auxiliary problems 3.2.1 Type I
As stated above, I a fields can be used as a basis forū and thus for u−u o . Then, using the condensation of the DIC problem (Section 2.2) as an auxiliary problem, the contribution of the I a fields is obtained. Once a parametrization of I a is defined (see Section 3.3), a condensed DIC problem is used to close the kinematic gapū.

Type II
Conversely, II a fields have no effect onū; they are not observable and thus cannot be obtained from DIC. To retrieve uniqueness, one has to formulate an auxiliary problem constraining II a . It is proposed in the following that some information are given concerning the thermodynamics of the deformation process. For illustration purposes, 3 cases are investigated: 1. The material is homogeneous, its behaviour is linear, and there is no dissipation: In this case, there exists  such that =  . Assuming that is known (from DIC), then the governing equation for II a is 2. The material is not homogeneous, but it is isotropic ( o is taken isotropic as well), its behaviour is linear, and there is no dissipation. In this case, one can define a contrast field such that =  o . Assuming that is known (from DIC), then the governing equation 3. The material is homogeneous, isotropic, and elastoplastic with isotropic hardening. Its elastic behaviour is defined by  o that is assumed to be isotropic. Then assuming that the plastic flow occurs following the normality rule, the principle of maximum dissipation applies. The problem is rephrased incrementally, and an increment Δ II a of anelastic strain between 2 time instants is searched to maximize Searching for the stationarity of the dissipation as defined above leads to a linear problem in Δ II a .
Obviously, these auxiliary problems cannot be solved pointwise. One first needs a parametrization of II a and then to solve these problems in a least-squares sense. The former points is the purpose of the next section.

Parametrization
Using usual FE settings, a 2D mesh of  e simplex elements is considered. From this mesh, a strain basis of 3 ×  e components is considered. The 3 ×  e FE internal force vectors F int (2 ×  components) with  o as constitutive operator are elaborated. A singular value decomposition of these vectors is performed to extract their linear combination giving rise to non-zero internal forces. Up to  I = 2 ×  − 3 combinations (denoted as E a ) can be obtained. If K o denotes the FE elastic stiffness matrix with  o as constitutive operator, then  I displacement solutions V are obtained from the resolution of the linear systems defined as The stress fields are computed, and type I and type II fields are computed according to their respective definitions (Equations 17 and (18)).
The following parametrizations are introduced: Note that the number of terms in each of these parametrizations can be adjusted independently, namely,  I and  II . For solving the auxiliary problem I (Section 3.2.1), one needs a parametrization of the displacement u. The first component of this parametrization is U o . This displacement vector is obtained from an FE elastic simulation using  o as the elastic tensor and a prescribed displacement measured by DIC along the boundary of the ROI where non-zero stress conditions apply. The amplitude of this displacement is adjusted so that the resulting forces F o match the experimental data. Note that because  o may not be the actual constitutive operator, after adjustment of F o , the displacement along the boundary of the ROI does not match the measured displacement anymore. The next terms of the displacement parametrization are obtained from the anelastic fields, namely, U I . This parametrization does not contain rigid body motions (2 translations and one rotation in 2D) that are considered through additional basis vectors in U R . These bases are then concatenated asŪ = [U I U R ] for describingū.

Resolution 3.4.1 Type I
A parametrization of the anelastic strain field is proposed above. Concerning the displacement, the amplitude of one of the basis functions is known a priori to ensure that the external force resulting from the measured stress matches the measured force. The deformed image g is thus advected by u o +ū and the DIC problem solved with respect to the parameter involved in the parametrization [Ū I ] only. Following Section 2.2, the DIC problem is solved with a condensation operator L equal toŪ I .
Using the proposed methodology, DIC provides for the  I amplitudes of [Ū I ] components, let us sayΛ I . This of course gives a direct access to the measured displacement and total strain field. In a post-processing step, anelastic strain field I a is derived from E I a and the corresponding components ofΛ I .

Type II
For computing II a , one of the auxiliary problems presented in Section 3.2.2 is solved. All these problems are linear, and once the parametrization E II a is adopted, they recast as where P is a  e ×  II matrix collecting the linear constraints of the auxiliary problem II in each element of the mesh and Q collects  e values for the constraints. The identity matrix I and the additional parameter's vector A account for the contribution of additional quantities like  for case 1 or unknowns defining the variation of for case 2. This linear system of equations is solved in a least-square sense providing for the  II amplitudesΛ II of the parametrization for II a .ī s then reconstructed using Σ II , and the actual stress field is reconstructed. Also, the additional parameters are identified like the actual elastic tensor  for case 1 or the contrast field for case 2.
Note that while thermodynamic assumptions considered in the auxiliary problems constrain the type of constitutive law satisfied by the actual fields, the decomposition adopted herein is completely generic. Another important point is that any combination of components of the proposed basis satisfies the local balance of momentum.

EXAMPLES
To illustrate the capability of the proposed method to measure stress, strain, and anelastic strain fields from digital images, synthetic data are generated. An FE simulation is run, and the result from this simulation is used to generate a set of deformed images given a reference speckle pattern. Three examples illustrating the 3 auxiliary problems formulated in the previous are proposed.

Anisotropic linear elasticity
In this example, the model consists in a rectangular domain of 1800 × 3800 pixels with 3 holes of different diameter. The reference image is shown by Figure 1A. In the FE simulation used for generating the synthetic data, the displacement along the bottom edge is clamped, while along the top edge, the vertical displacement is set to 0, 30, and 30 pixels. The horizontal displacement is prescribed to 30, 0, and 30 pixels. In this numerical simulation, the material is assigned a linear orthotropic behaviour  i defined as in the x, y coordinate system of the reference image. The synthetic images are generated using a very fine mesh, 4 times finer than that used in the following for MIC, and depicted in Figure 1B. In this mesh, the element size is 60 pixels far For the MIC analysis,  o is set as an isotropic elastic tensor having its Young modulus and Poisson ratio set to 5 GPa and 0.2, respectively. The number of I a functions is set as high as possible (2 × 1539 − 3) whereas II a is described with 500 modes. The displacement fields for the 3 load steps are presented in Figure 2. For illustration purposes, the reference stress field σ o and the reconstructed stress field σ are plotted in Figure 3 for the last load step. It is worth noting that they differ in the vicinity of the holes while they are similar close to the bottom and top edges as their vertical components are the main contribution to the external load that both of them balance. To assess the reliability of the stress estimate, it is compared on Figure 4 to a stress field obtained by applying the actual elastic tensor  i to the measured strain field. Indeed, due to discretization effects, it is not possible to compare the reconstructed stress to the numerical simulation used for image generation. The error in the stress field is plotted in Figure 4 with different colour scales. On Figure 4B, the colour scale has the same amplitude as for representing σ itself whereas this amplitude is reduced by a factor of 10 in Figure 4A. This allows one to appreciate the robustness of the estimation. Indeed, the average difference between the estimated stress field and its theoretical value is about 9 Mpa and the standard deviation of this difference is about 67 Mpa. Compared to a maximum stress value higher than 1000 Mpa, it is stated that the uncertainty is around 5%. Also, the auxiliary problem provides an estimate of the elastic tensor. The following values are obtained: GPa. (28) The accuracy is again around 5%. To further illustrate the capability of the proposed approach, a 3-dimensional plot of the vertical stress as a function of the measured vertical and horizontal strain is proposed in Figure 5. One might notice that the region of the stress/strain domain with low signal (low stress and low vertical strain) concentrates the points with highest error. There is also high error points for high stress corresponding to the hole boundary where the strain measurement is less reliable. However, it is shown through this example that MIC provides not only for displacement and strain fields but also for reliable stress field estimates as well as the full elastic tensor from one single test configuration.

Heterogeneous isotropic elasticity
In this second example, the same rectangular plate as before is considered. However, instead of having holes, a soft circular inclusion (the contrast is set to 0.1) is included in the centre of the plate (see Figure 6A). The same loading scenario as in the previous example is modelled. While the mesh used for the simulation generating the images explicitly describes the inclusion, the mesh for MIC and the images are built without considering the inclusion (see Figure 6B and 6C). The MIC mesh contains 1969 nodes for 3773 elements of around 50 pixels. For the MIC analysis,  o is set as an isotropic elastic tensor having its Young modulus and Poisson ratio set to 5 GPa and 0.2, respectively. This is the elastic model used outside the inclusion for generating the images. The number of modes for I a is set to 500, while for II a , it is set to 100 modes. The variation of the contrast field is described with the FE linear shape function of the mesh in Figure 6C. It is worth noting that the material discontinuity considered for generating the images cannot be captured exactly. A smoothing of the identified contrast field is expected across the inclusion boundary.
On Figure 7, the vertical component of the strain tensor for the 3 load steps are presented. One clearly observes the strain concentration in the matrix around the hole and very high but quasi-homogeneous strain within the soft inclusion. Using the auxiliary problem 2 in Section 3.2.2, II a is obtained for reconstructing the stress field from the reference solution. As the material is homogeneous in the reference problem, for the second load step that corresponds to uniaxial tension, σ o is almost homogeneous in the sample (see Figure 8A). The gap to be compensated bȳto recover the actual stress is large and heterogeneous. However, in Figure 8B, the expected stress concentration around the hole is recovered. The emphasis is here on the ability of the approach to identify material property field. For this purpose, the identified contrast field is plotted on Figure 8C. The white circle in this figure defines the actual boundary of the inclusion as considered in the simulation used for generating the images. Although the mesh used to perform MIC is uniform, a very good agreement is obtained on the geometry of the inclusion. Further, the stiffness contrast between the inclusion and the matrix is captured. The average error on the contrast over the whole sample is 0.02, while the standard deviation of the error is 0.13.

Elastoplasticity
A rectangular sample of 100×20 mm 2 with a hole of diameter 5 mm is considered. It is assigned an elastoplastic behaviour with linear isotropic hardening. The Young modulus is 168 GPa, the initial yield stress 284 MPa, and the hardening coefficient 1.48 GPa. The Poisson ratio is 0.25. The numerical simulation is run with a mesh similar to that shown in Figure 9A but with one additional subdivision. A 2D plane stress formulation is used. The reference image is depicted in Figure 9B and last deformed image in Figure 9C. The size of the images is 5120 × 1400 pixels, corresponding to a physical pixel size of 22.2 m. The speckle pattern for the initial image is a real one taken from an image of a real sample. The load history applied in the simulation is illustrated in Figure 10. Two hundred steps are performed. Two unloading stages are considered at steps 40 and 120. The MIC analysis is run using 500 type I modes and 300 type II modes. The elastic tensor  o is assigned the same value as the elastic tensor for the numerical simulation used to generate the images. In the DIC analysis, the displacement estimation along the boundary of the FE mesh is usually less robust than in the bulk. To avoid that this lower robustness affects the result of the MIC analysis where the DIC displacement is used as boundary conditions for computing the admissible elastic solution u o , one layer of elements is removed at both ends of the mesh. The DIC displacement used    Figure 9A for the last loading step are depicted in Figure 11. The scatter between the 2 displacements is difficult to observe by the naked eyes. For assessing the accuracy of the displacement field obtained by MIC, the correlation error map is plotted for the 2 analyses in Figure 12. These 2 error maps are again difficult to distinguish, what confirms that the MIC displacement leads to a registration of the images of the same accuracy as the DIC analysis. To further assess this point, the average correlation error for DIC and MIC is plotted as a function of the steps in Figure 13. The error level is slightly higher for MIC, what is expected because only 500 modes are used to measure the displacement instead of 2 ×  . The search of the optimal solution in terms of grey-level conservation is thus more constraint, and the minimum error level reached with the reduced parametrization is higher. These results ensure that the displacement parametrization adopted for MIC leads to measured fields globally as accurate as those obtained for a DIC analysis.
Mechanical image correlation also provides for the estimation of meaningful mechanical fields. First, the evolution of the stress component in the loading direction is provided in Figure 14. The decrease in stress during unloading is well captured, as well as its progressive increase during the loading stages.
In Figure 15, the norm of the anelastic strain̄a is plotted for the same steps as the stress in Figure 14. It is observed, as expected for an elastoplastic material for which a corresponds to the plastic strain, that the anelastic strain is kept  To further illustrate the capability of MIC to measure stress fields, an analysis of the entire set of measured material states is provided in Figure 17. This set consists in more than 600 000 points in the space defined by the components of stress, strain, and anelastic strain tensors. For visualization purposes, one needs to reduce the dimensionality of this database. Because it is assumed that von Mises plasticity holds, the von Mises stress is computed as well as the deviatoric part of the elastic strain and the norm of the anelastic strain. In Figure 17A, the evolution of the von Mises stress is plotted as a function of the norm of the plastic strain for a few points. The typical loading unloading scenario of elastoplasticity is captured. In Figure 17B, von Mises stress is plotted as a function of the deviatoric part of the elastic strain and the norm

CONCLUSION
A methodology is proposed to build, in a general case, a parametric description of the displacement and stress fields based on an anelastic strain field. Digital image correlation is used to obtained the displacement parameters whereas the stress parameters are obtained thanks to thermodynamic assumptions on the deformation process leading to an auxiliary problem. The method thus allows for measuring displacement, strain, and stress fields. Further, an anelastic strain field measuring the distance between the actual local constitutive behaviour and elasticity according to a given but arbitrary elastic tensor is obtained. The ability of the method to measure stress fields in different situations is illustrated. The raw data extracted from the measurement are admissible material states in a space of dimension 9 (3 components for the strain, stress, and anelastic strain tensors) for 2D problems.