# Computational measurements of stress fields from digital images

Julien Réthoré, Adrien Leygue, Michel Coret, Laurent Stainier, Erwan Verron

► To cite this version:

Julien Réthoré, Adrien Leygue, Michel Coret, Laurent Stainier, Erwan Verron. Computational measurements of stress fields from digital images. International Journal for Numerical Methods in Engineering, 2018, 113 (12), pp.1810-1826. 10.1002/nme.5721. hal-01454432v3

HAL Id: hal-01454432

https://hal.science/hal-01454432v3

Submitted on 12 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Computational measurements of stress fields from digital images

Julien Réthoré | Adrien Leygue | Michel Coret | Laurent Stainier | Erwan Verron

Research Institute in Civil and Mechanical Engineering (GeM), Ecole Centrale de Nantes, Université de Nantes, UMR 6183 CNRS, 1 rue de la Noë, F-44 321 Nantes, France

**Correspondence**

Julien Réthoré, Research Institute in Civil and Mechanical Engineering (GeM), Ecole Centrale de Nantes, Université de Nantes, UMR 6183 CNRS, 1 rue de la Noë, F-44 321 Nantes, France.
Email: Julien.Rethore@ec-nantes.fr

**Funding information**

Région Pays de la Loire and Nantes Métropole

**Summary**

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.

**KEYWORDS**

digital image correlation, identification, stress measurements

## 1 | INTRODUCTION

Since the pioneering work of Lucas, Kanade, and others,¹ 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²) because they are the techniques implemented in commercial softwares. One should also mention the non-parametric formulation proposed by Vercauteren et al³ 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.⁴ 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⁵ 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.⁶,⁷ 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⁸ 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⁹ 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.

# 2 | DIGITAL IMAGE CORRELATION

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

## 2.1 | 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:

$$ \mathbf{u} = \operatorname{Arg Min} \sum_{p \in \text{ROI}} [f(\mathbf{x}_p) - g(\mathbf{x}_p + \mathbf{u}(\mathbf{x}_p))]^2, \quad (1) $$

where *p* is the index of a pixel within the region of interest (ROI) over which the computation is performed and **x**<sub>p</sub> 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.³ 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:

$$ \mathbf{u}(\mathbf{x}) = \sum_{k \in \mathcal{N}} u_k \mathbf{N}_k(\mathbf{x}). \quad (2) $$

In this decomposition, *N* is the set of degrees of freedom (DOFs) *u*<sub>*k*</sub>, and *N*<sub>*k*</sub> 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.¹⁰
Other specific bases, like for 4-point bending test¹¹ or Brazilian test,¹² can be used, but for arbitrary displacement, FE-based⁴,¹³ or CAD (Computer-Aided Design)-based kinematics¹⁴,¹⁵ 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

$$[\tilde{M}] \{dU\} = \{\tilde{b}\}, \qquad (3)$$

with

$$[\tilde{M}] = [N]^T [\nabla G] [\nabla G] [N], \qquad (4)$$

and

$$\{\tilde{b}\} = [N]^T [\nabla G] \{F - G\}, \qquad (5)$$

$\nabla$ denoting spatial derivation. [N] is a matrix that collects the value of the basis functions at the pixels, [$\nabla G$] a diagonal matrix that collects the values of $\nabla 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.¹⁶ To reduce computational cost, a modified algorithm is adopted: We substitute $\nabla G$, which depends on the current solution vector, by $\nabla F$. The modified system to solve thus becomes

$$[M] \{dU\} = \{b\}, \qquad (6)$$

with

$$[M] = [N]^T [\nabla F] [\nabla F] [N], \qquad (7)$$

and

$$\{b\} = [N]^T [\nabla F] \{F - G\}. \qquad (8)$$

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, $\nabla F$ equals $\nabla G$. After a few iterations (2 or 3), the modified algorithm and the initial one reveal close to identical.

## 2.2 | 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,¹⁰ a direct access to stress intensity factors is provided. An other interesting example is described in Leclerc et al⁶ and Réthoré⁷: 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

$$\{U\} = [L] \{\Lambda\}. \qquad (9)$$

In this formalism, $\Lambda$ 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:

$$[M_C] \{d\Lambda\} = \{b_C\}, \qquad (10)$$

with

$$[M_C] = [L]^T [M] [L], \quad \{b_C\} = [L]^T \{b\}. \qquad (11)$$
3 | MECHANICAL IMAGE CORRELATION

The approaches in Leclerc et al⁶ and Réthoré⁷ are parametric identification techniques: They allow for estimating experi-
mentally 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 evi-
denced 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.

**3.1 | Formulation**

Let $C_o$ be an elastic operator. Then, an elastic solution $\mathbf{u}_o, \sigma_o$ according to $C_o$ and balancing the measured load $F_o$ is computed using a standard FE solver. The vector collecting the nodal displacements of this elementary solution is denoted by $\mathbf{U}_o$. Because $C_o$ differs from the actual constitutive law of the material, this reference solution $\mathbf{u}_o, \sigma_o$ does not match the actual one $\mathbf{u}, \sigma$. One can define a gap in terms of displacement (kinematic gap) $\bar{\mathbf{u}}$ such that

$$
\mathbf{u} = \mathbf{u}_o + \bar{\mathbf{u}} \tag{12}
$$

and a gap in terms of stress (static gap) $\bar{\sigma}$

$$
\sigma = \sigma_o + \bar{\sigma}. \tag{13}
$$

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 $\sigma_o$ and $\sigma$ balance the external load, $\bar{\sigma}$ is a self-balanced stress field. It admits the following properties:

$$
\operatorname{div} \bar{\sigma}=0 \quad \text { over the ROI, }
$$

However, there is a priori no straightforward relation between $\bar{\sigma}$ and $\bar{\mathbf{u}}$. Without loss of generality, one introduces an
anelastic strain field $\epsilon_a$ such that

$$
\bar{\sigma} = C_0(\bar{\epsilon} - \epsilon_a), \quad (15)
$$

where $\bar{\epsilon}$ is the symmetric gradient of the kinematic gap $\bar{\mathbf{u}}$. Then using Equation 14, one obtains

$$
\operatorname{div}(C_o \bar{\epsilon}) - \operatorname{div}(C_o \epsilon_a) = 0 \quad \text{over the ROI,} \quad (C_o(\bar{\epsilon} - \epsilon_a)n = 0 \quad \text{over the boundary of ROI.} \quad (16)
$$

It is thus clear that if $\epsilon_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
$\epsilon_a$, and second, find a parametrization of $\epsilon_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 $\epsilon_a$ fields in 2 contributions:

$$
\epsilon_a = \epsilon_a^I + \epsilon_a^{II}. \tag{17}
$$

Once the problem defined by Equation 16 is solved, a solution $\bar{\mathbf{u}}$ is obtained up to a rigid body motion. It is then possible
to define the compatible part $\epsilon_a^I$ of $\epsilon_a$ as the symmetric gradient of $\bar{\mathbf{u}}$:

$$
\epsilon_a^I = \bar{\epsilon} \tag{18}
$$

The remaining part $\epsilon_a^{II} = \epsilon_a - \epsilon_a^I$ is not compatible in general, and it is not associated to any displacement ($\bar{\mathbf{u}}^{II} = 0$) as

$$
\operatorname{div}(C_o \epsilon_a^{II}) = 0. \tag{19}
$$

The resulting self-balanced stress field is then written as

$$
\bar{\sigma}^{II} = -C_o \epsilon_a^{II}. \tag{20}
$$

Conversely, $\epsilon_a^I$ induces a displacement $\bar{\mathbf{u}}_a^I$ but no stress as $\epsilon_a^I = \bar{\epsilon}$.

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:

* $\epsilon_a^I$ being compatible, it is intrinsically related to a displacement (namely, $\bar{\mathbf{u}}$). If a set of basis functions deriving from a set displacement functions is chosen for $\epsilon_a^I$, they can be used as a basis for $\mathbf{u}$ by using Equation 12.

* the effect of $\epsilon_a^{II}$ is not observable ($\bar{\mathbf{u}}^{II} = 0$).
## 3.2 | Auxiliary problems

### 3.2.1 | Type I

As stated above, $\epsilon_a^I$ fields can be used as a basis for $\bar{u}$ 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 $\epsilon_a^I$ fields is obtained. Once a parametrization of $\epsilon_a^I$ is defined (see Section 3.3), a condensed DIC problem is used to close the kinematic gap $\bar{u}$.

### 3.2.2 | Type II

Conversely, $\epsilon_a^{II}$ fields have no effect on $\bar{u}$; they are not observable and thus cannot be obtained from DIC. To retrieve uniqueness, one has to formulate an auxiliary problem constraining $\epsilon_a^{II}$. 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 $C$ such that $\sigma = C\epsilon$. Assuming that $\epsilon$ is known (from DIC), then the governing equation for $\epsilon_a^{II}$ is

$$C_o\epsilon_a^{II} = \sigma_o - C\epsilon. \qquad (21)$$

2. The material is not homogeneous, but it is isotropic ($C_o$ is taken isotropic as well), its behaviour is linear, and there is no dissipation. In this case, one can define a contrast field $\phi$ such that $\sigma = \phi C_o \epsilon$. Assuming that $\epsilon$ is known (from DIC), then the governing equation for $\epsilon_a^{II}$ is

$$C_o\epsilon_a^{II} = \sigma_o - \phi C_o \epsilon. \qquad (22)$$

3. The material is homogeneous, isotropic, and elastoplastic with isotropic hardening. Its elastic behaviour is defined by $C_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 $\Delta\epsilon_a^{II}$ of anelastic strain between 2 time instants is searched to maximize

$$(\sigma_o - C_o(\epsilon_a^{II} + \Delta\epsilon_a^{II})) : (\Delta\epsilon_a^I + \Delta\epsilon_a^{II}). \qquad (23)$$

Searching for the stationarity of the dissipation as defined above leads to a linear problem in $\Delta\epsilon_a^{II}$.

Obviously, these auxiliary problems cannot be solved pointwise. One first needs a parametrization of $\epsilon_a^{II}$ and then to solve these problems in a least-squares sense. The former points is the purpose of the next section.

## 3.3 | Parametrization

Using usual FE settings, a 2D mesh of $N_e$ simplex elements is considered. From this mesh, a strain basis of $3 \times N_e$ components is considered. The $3 \times N_e$ FE internal force vectors $\mathbf{F}_{int}$ ($2 \times N$ components) with $C_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 $N_I = 2 \times N - 3$ combinations (denoted as $\mathbf{E}_a$) can be obtained. If $\mathbf{K}_o$ denotes the FE elastic stiffness matrix with $C_o$ as constitutive operator, then $N_I$ displacement solutions $\mathbf{V}$ are obtained from the resolution of the linear systems defined as

$$[\mathbf{K}_o] \{\mathbf{V}\} = [\mathbf{F}_{int}]\{\mathbf{E}_a\}. \qquad (24)$$

The stress fields

$$\Sigma = C_o(\nabla^S \mathbf{V} - \mathbf{E}_a) \qquad (25)$$

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:

$$\text{I} \quad \bar{\mathbf{U}}^I = \mathbf{V}, \mathbf{E}_a^I = \nabla^S \mathbf{V}, \Sigma^I = 0,$$

$$\text{II} \quad \mathbf{E}_a^{II} = \mathbf{E}_a - \nabla^S \mathbf{V}, \bar{\mathbf{U}}^{II} = 0, \Sigma^{II} = \Sigma.$$

Note that the number of terms in each of these parametrizations can be adjusted independently, namely, $N_I$ and $N_{II}$. For solving the auxiliary problem I (Section 3.2.1), one needs a parametrization of the displacement $\mathbf{u}$. The first
component of this parametrization is $\mathbf{U}_o$. This displacement vector is obtained from an FE elastic simulation using $C_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 $C_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, $\mathbf{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 $\mathbf{U}_R$. These bases are then concatenated as $\bar{\mathbf{U}} = [\mathbf{U}^I \mathbf{U}_R]$ for describing $\bar{\mathbf{u}}$.

## 3.4 | 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 $\mathbf{u}_o+\bar{\mathbf{u}}$ and the DIC problem solved with respect to the parameter involved in the parametrization $[\bar{\mathbf{U}}^I]$ only. Following Section 2.2, the DIC problem is solved with a condensation operator $\mathbf{L}$ equal to $\bar{\mathbf{U}}^I$.

Using the proposed methodology, DIC provides for the $\mathcal{N}_I$ amplitudes of $[\bar{\mathbf{U}}^I]$ components, let us say $\bar{\mathbf{\Lambda}}^I$. This of course gives a direct access to the measured displacement and total strain field. In a post-processing step, anelastic strain field $\epsilon_a^I$ is derived from $\mathbf{E}_a^I$ and the corresponding components of $\bar{\mathbf{\Lambda}}^I$.

### 3.4.2 | Type II

For computing $\epsilon_a^{II}$, one of the auxiliary problems presented in Section 3.2.2 is solved. All these problems are linear, and once the parametrization $\mathbf{E}_a^{II}$ is adopted, they recast as

$$[P][\mathbf{E}_a^{II} \mathbf{I}]{\bar{\mathbf{\Lambda}}^{II} \mathbf{A}} = \{Q\}, \qquad (26)$$

where $P$ is a $\mathcal{N}_e \times \mathcal{N}_{II}$ matrix collecting the linear constraints of the auxiliary problem $II$ in each element of the mesh and $Q$ collects $\mathcal{N}_e$ values for the constraints. The identity matrix $\mathbf{I}$ and the additional parameter's vector $\mathbf{A}$ account for the contribution of additional quantities like $C$ for case 1 or unknowns defining the variation of $\phi$ for case 2. This linear system of equations is solved in a least-square sense providing for the $\mathcal{N}_{II}$ amplitudes $\bar{\mathbf{\Lambda}}^{II}$ of the parametrization for $\epsilon_a^{II}$. $\bar{\mathbf{\sigma}}$ is then reconstructed using $\Sigma^{II}$, and the actual stress field is reconstructed. Also, the additional parameters are identified like the actual elastic tensor $C$ for case 1 or the contrast field $\phi$ 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.

# 4 | 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.

## 4.1 | 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 Cᵢ defined as

$$C_i = \begin{bmatrix} 5.56 & 1.39 & 0.00 \\ 1.39 & 55.56 & 0.00 \\ 0.00 & 0.00 & 4.16 \end{bmatrix}_{x,y} \text{ GPa,} \qquad (27)$$

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

<figure><img src="image_2.png" /><figcaption>FIGURE 1 Reference image (1800 × 3800 pixels) and mesh for mechanical image correlation analysis in the case of anisotropic elasticity</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2 Displacement magnitude in pixel for the 3 load steps in the case of anisotropic elasticity (amplification factor 5)</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**FIGURE 3** Reference stress and stress field estimate (vertical component) in the case of anisotropic elasticity (last load step, amplification factor 5)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**FIGURE 4** Error assessment on the vertical component of the estimated stress in the case of anisotropic elasticity (last load step, amplification factor 5). A, With the same colour scale as for Figure 3. B, With a lower data range in the colour scale</figcaption></figure>

from the holes and 40 pixels in the vicinity of the holes. This mesh contains 1539 nodes and 2808 triangular elements with linear interpolation.

For the MIC analysis, $C_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 $\epsilon_a^l$ functions is set as high as possible (2 × 1539 − 3) whereas $\epsilon_a^{ll}$ 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 $\sigma_o$ and the reconstructed stress field $\sigma$ 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 $C_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:

$$ C = \begin{bmatrix} 7.4 & 1.60 & 0.05 \\ 1.60 & 51.75 & 0.65 \\ 0.05 & 0.65 & 4.60 \end{bmatrix}_{xy} \text{ GPa.} \qquad (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.

<figure><img src="image_6.png" /><figcaption>FIGURE 5 Vertical stress as a function of the measured vertical and horizontal strain in the case of anisotropic elasticity. The dots correspond to measurements, while the plane is the prescribed elastic tensor. The colour of the dots defines the distance to the plane</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6 Test configuration, reference image (1800 × 3800 pixels), and mesh for mechanical image correlation analysis in the case of heterogeneous elasticity</figcaption></figure>
## 4.2 | 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, $C_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 $\epsilon_a^I$ is set to 500, while for $\epsilon_a^{II}$, it is set to 100 modes. The variation of the contrast field $\phi$ 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, $\epsilon_a^{II}$ 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, $\sigma_o$ is almost homogeneous in the sample (see Figure 8A). The gap to be compensated by $\bar{\sigma}$ to recover the actual stress is large and heterogeneous. However, in Figure 8B, the expected stress concentration around the hole is recovered.

<figure><img src="image_8.png" /><figcaption>FIGURE 7 Vertical strain field for the 3 load steps in the case of heterogeneous elasticity (amplification factor 5)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 8 Reference stress and stress field estimate (vertical component) and identified contrast field (second step, amplification factor 5)</figcaption></figure>
The emphasis is here on the ability of the approach to identify material property field. For this purpose, the identified contrast field $\phi$ 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.

## 4.3 | Elastoplasticity

A rectangular sample of 100×20 mm² 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 $C_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><img src="image_10.png" /><figcaption>**FIGURE 9** Finite element mesh for the mechanical image correlation analysis. Initial and last images generated from the finite element simulation (5120 × 1400 pixels)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**FIGURE 10** Global loading history for the finite element simulation</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 11** Displacement (in pixel) in the loading direction at the last step. A, For a digital image correlation analysis. B, For a mechanical image correlation analysis</figcaption></figure>

as boundary conditions is thus picked up from nodes in the bulk of the DIC mesh. The displacement fields measured by DIC and MIC using the mesh in 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 \times N$. The search of the optimal solution in terms of grey-level conservation is thus more
<figure><img src="image_13.png" /><figcaption>**FIGURE 12** Correlation error map in grey level at the last step. A, For a digital image correlation analysis. B, For a mechanical image correlation analysis</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 13** Comparison of the average correlation error in grey level for the digital image correlation (DIC) and mechanical image correlation (MIC) analysis</figcaption></figure>

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 $\epsilon_a$ is plotted for the same steps as the stress in Figure 14. It is observed, as expected for an elastoplastic material for which $\epsilon_a$ corresponds to the plastic strain, that the anelastic strain is kept
<figure><img src="image_15.png" /><figcaption>**FIGURE 14** Evolution of the stress component in the axial direction in MPa. A, Just before the first unloading. B, At the first unloading. C, Just before the second unloading. D, At the second unloading. E, At the last step</figcaption></figure>

constant during unloading. Note that the norm of the tensor includes the contribution of the out-of-plane component that has been derived using plastic incompressibility. The latter is plotted for the last step in Figure 16. Expected results like shrinkage at the pole of the hole are retrieved.
<figure><img src="image_16.png" /><figcaption>**FIGURE 15** Evolution of the anelastic strain norm. A, Just before the first unloading. B, At the first unloading. C, Just before the second unloading. D, At the second unloading. E, At the last step</figcaption></figure>

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
<figure><img src="image_17.png" /><figcaption>**FIGURE 16** Out-of-plane component of the plastic strain at the last step</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**FIGURE 17** Stress-strain plots for the material state database obtained by mechanical image correlation</figcaption></figure>

of the plastic strain. The shape of the domain of admissible material states agrees well with the one considered in the
numerical simulations to generate the images.

5 | 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 obtain 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.

ACKNOWLEDGEMENT

The support of Région Pays de la Loire and Nantes Métropole through grant Connect Talent IDS is gratefully acknowl-
edged.

ORCID

Julien Réthore ID http://orcid.org/0000-0003-0690-8549
Laurent Stainier ID http://orcid.org/0000-0001-6719-6616

REFERENCES

1. Lucas BD, Kanade T. An iterative image registration technique with an application to stereo vision. *IJCAI*. 1981;81:674-679.
2. Sutton M, Cheng M, Peters W, Chao Y, McNeill S. Application of an optimized digital image correlation method to planar deformation analysis. Image Vision Comput. 1986;4(3):143-150.

3. Vercauteren T, Pennec X, Perchant A, Ayache N. Diffeomorphic demons: efficient non-parametric image registration. NeuroImage. 2009;45(1):S61-S72.

4. Besnard G, Hild F, Roux S. 'Finite-element' displacement fields analysis from digital images: application to Portevin-Le Châtelier bands. Exp Mech. 2006;46(6):789-803.

5. Réthoré J, Hild F, Roux S. An extended and integrated digital image correlation technique applied to the analysis of fractured samples. Eur J Comput Mech. 2009;18:285-306.

6. Leclerc H, Perie J, Roux S, Hild F. Computer Vision/Computer Graphics Collaboration Techniques, chap. Integrated Digital Image Correlation for the Identification of Mechanical Properties. Berlin: Springer; 2009.

7. Réthoré J. A fully integrated noise robust strategy for the identification of constitutive laws from digital images. Int J Numer Methods Eng. 2010;84:631-660. https://doi.org/10.1002/nme.2908.

8. van der Maaten L, Postam E, van den Herik J. Dimensionality reduction: a comparative review. Technical Report, Tilburg, Netherlands: Tilburg Centre for Creative Computing, Tilburg University, Technical Report: 2009-005; 2009.

9. Kirchdoerfer T, Ortiz M. Data-driven computational mechanics. Comput Methods Appl Mech Eng. 2016;304:81-101.

10. Roux S, Hild F. Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches. Int J Fract. 2006;140(1-4):141-157.

11. Leplay P, Réthoré J, Meille S, Baietto M. Damage law identification of a quasi brittle ceramic from a bending test using digital image correlation. J Eur Ceram Soc. 2010;30(13):2715-2725.

12. Hild F, Roux S. Digital image correlation: from displacement measurement to identification of elastic properties—a review. Strain. 2006;42(2):69-80.

13. Réthoré J, Hild F, Roux S. Extended digital image correlation with crack shape optimization. Int J Numer Methods Eng. 2008;73(2):248-272. https://doi.org/10.1002/nme.2070.

14. Elguedj T, Réthoré J, Buteri A. Isogeometric analysis for strain field measurements. Comput Methods Appl Mech Eng. 2011;200(1-4):40-56.

15. Réthoré J, Elguedj T, Simon P, Coret M. On the use of nurbs functions for displacement derivatives measurement by digital image correlation. Exp Mech. 2009;50(7):1099-1116. https://doi.org/10.1007/s11340-009-9304-z.

16. Réthoré J, Hild F, Roux S. Shear-band capturing using a multiscale extended digital image correlation technique. Comput Methods Appl Mech Eng. 2007;196(49-52):5016-5030. https://doi.org/10.1016/j.cma.2007.06.019.