# Inverse identification of the hardening behavior with a fully 3D virtual fields method

Miguel G Oliveira, Sandrine Thuillier, António Andrade-Campos

► To cite this version:

Miguel G Oliveira, Sandrine Thuillier, António Andrade-Campos. Inverse identification of the hardening behavior with a fully 3D virtual fields method. Advanced Modeling and Simulation in Engineering Sciences, 2025, 12 (1), pp.10. 10.1186/s40323-025-00293-7. http://05050227

HAL Id: hal-05050227

https://hal.science/hal-05050227v1

Submitted on 29 Apr 2025

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>
RESEARCH ARTICLE

Open Access

<figure><img src="image_2.png" /></figure>

Check for
updates

# Inverse identification of the hardening behavior with a fully 3D virtual fields method

Miguel G. Oliveira<fnref n="1" />,<fnref n="2" />*, Sandrine Thuillier<fnref n="2" /> and António Andrade-Campos<fnref n="1" /> ID

## Abstract

This work investigates the use of the virtual fields method (VFM) to calibrate the hardening behavior of sheet metals up to rupture. Indeed, a few studies dealing with the VFM in a 2D framework have shown that after necking, the validity of the VFM is compromised, leading to inaccurate descriptions of the material hardening behavior. However, no study has yet quantified the validity of the VFM in a 2D framework. Therefore, it is essential to assess its limitations and, in turn, highlight the advantages of a 3D approach. This study investigates the influence of sample thickness on both 2D and 3D VFM, first by analyzing the evolution of the predicted load and then by evaluating the identification of material parameters in the Swift hardening law. In the first part, the 3D-VFM has shown to yield accurate results while using the 2D-VFM resulted in non-negligible relative errors after necking, showing that the error is correlated to the thickness. In the second part, the identification of the parameters of the Swift law was carried out using different combinations of three manually defined virtual fields. The material parameters identified with the 2D-VFM did not match the reference parameters while the 3D-VFM was able to correctly identify them in every instance.

**Keywords:** Virtual fields method, Parameter identification, Three-dimensional framework, Hardening behavior, Necking

## Introduction

Conventionally, the elastoplastic behavior of sheet metals is characterized through standard uniaxial tensile tests. From these tests, it is possible to extract the true stress–strain curve essential to characterize the strain hardening behavior. The former is important because many sheet metal forming processes achieve plastic strains beyond the point of maximum uniform elongation. Moreover, to simulate these processes in finite element analysis (FEA) software, the key requirement is the definition of the true stress–strain curve. Without an accurate characterization of this curve, numerical simulations will not be representative of the real material behavior.

However, diffuse necking is observed after the point of maximum uniform elongation in the standard uniaxial tensile tests [1–3]. The diffuse necking leads to a non-uniform

<figure><img src="image_3.png" /></figure>

<footnote>*Correspondence:
oliveiramiguel@ua.pt</footnote>

<footnote n="1">¹Centre for Mechanical
Technology and Automation
(TEMA), Department of
Mechanical Engineering,
University of Aveiro, 3810-193
Aveiro, Portugal</footnote>

<footnote n="2">²Univ. Bretagne Sud, UMR CNRS
6027, IRDL, F-56100 Lorient,
France</footnote>
deformation concentrated within the necking region [4,5], which does not allow measur-
ing the deformation through a conventional extensometer. Commonly, localized necking
is also observed in sheet metals after diffuse necking. While the former phenomenon
is associated with a decrease in the specimen width, in the latter, a rapid decrease of the
thickness is observed. A triaxial stress state also develops within this necking region, which
eventually leads to fracture of the material [6,7].

By measuring deformation with an extensometer, it will underestimate the resulting
strain values after the onset of necking. Therefore, several studies have proposed solutions
to determine the true stress–strain curve beyond this point [8]. Usually, the hardening
behavior is extrapolated from the pre-necking stress–strain curve, but can lead to large
errors. Bridgman [9] was one of the first authors to tackle the problem of characterizing
the post-necking hardening behavior. In this pioneering study, the author proposed an
analytical method based on geometrical parameters of the necking profile to correct and
compensate for necking in a round bar. Although this method has shown to perform well
for round tensile specimens, it requires significant experimental work and could not be
applied to specimens of rectangular cross-sections. However, rounded specimens are not
practical for thin sheets applications. Later, Zhang et al. [10] proposed a method that can be
applied to specimens of rectangular cross-section by determining a relationship between
the area and thickness reduction of the cross-section. A drawback of this method is the
required use of small specimens with thin sheets. Also, Mirone [11] proposed a material-
independent method that characterizes the necking behavior with reduced errors and
much less experimental effort compared to the method of Bridgman.

With the advance of full-field measurement techniques combined with inverse method-
ologies, it is possible the development of complete solutions that consider in the identifi-
cation procedure, the stress and strain at all measured material points lying in the surface
of the material. Most proposed solutions use a finite element (FE) based inverse method
for the identification of the post-hardening behavior [6,7,12–14]. Nevertheless, building
a FE model that accurately predicts the occurrence of necking is a hard task, for which
several aspects need to be considered. The FE models should be capable of predicting
the location of diffuse necking, which is often also a difficulty in real experiments. Special
attention should also be paid to building a FE mesh to avoid highly distorted elements, due
to the large strains involved in post-necking [15]. Employing information about the test’s
actual boundary conditions in the FE model can lead to more accurate results [16,17], but
getting this information is not straightforward. Another drawback of such methodologies,
which is often pointed out, is their high computational cost [18].

More recently, alternatives to FE based methodologies are becoming more popular by
avoiding their shortcomings and directly using the measured strain fields. Rossi et al.
[19] proposed an identification method, known as linear stress–strain curve identifica-
tion method (LSSCI), that uses a piecewise linear function to describe the stress–strain
curve even beyond necking. This method has been successfully applied, provided that high
quality full-field measurements are obtained up to the free edges of the specimen. Never-
theless, to fully characterize the elastoplastic behavior (strain hardening and anisotropy)
of sheet metals, a decoupled identification approach is required. In this approach, the
LSSCI method is used to characterize the hardening behavior, and the VFM is used to
identify the parameters of a yield criterion [20]. Coppieters et al. [21] proposed an energy
method, based on the observation that under quasi-static conditions, the internal work is
equal to the external work, assuming the hardening behavior is correctly described. Using
real experiments, it has been shown that this method is computationally efficient, without
the need to perform any FE simulations [22].

The LSSCI method is not fundamentally different from the VFM and can be considered
a special case of it. The VFM is based on the principle of virtual work and employs virtual
displacement fields instead of actual ones to balance the internal virtual work with the
external virtual work [23]. It is one of the most prominent inverse methodologies and
has already been applied to various applications [24–30]. In particular, it has already
been successfully applied in the characterization of the strain hardening behavior of thin
sheets up to large deformations [31]. However, some studies have highlighted that the
observed deviation between internal and external virtual works after diffuse necking is
primarily due to the 2D assumption [32], violating the validity of the VFM, leading to
an inaccurate description of the hardening behavior. The triaxial stress state can explain
the observed deviation that develops in the necking region [33]. This rich information
about the necking behavior is not used in the VFM because of the plane stress limitations
imposed by the lack of available three-dimensional (3D) full-field measurements. In fact,
in sheet metal characterization, the VFM is limited to be used in combination with surface
measurements, as techniques are not yet available to capture the 3D behavior inside the
material. To circumvent this limitation, the internal mesh generation method investigated
proposed by Rossi et al. [34] and investigated by Oliveira et al. [35] might provide the
necessary 3D displacement fields for using the VFM in a 3D framework. However, this
topic has not been extensively investigated, and at the best of the author’s knowledge, only
one study is found that extend the use of the VFM to a 3D framework [36], assuming no
plane stress or plane strain conditions.

Considering the amount of studies already using the VFM in the identification of mate-
rial parameters, its interest and potential is undeniable. Yet, it lacks a better understand-
ing of the consequences of its limitations, particularly the use of a two-dimensional (2D)
instead of a 3D framework. Although several problems might exist in the latter, such as the
lack of volumetric full-field measurements, it might encourage the scientific community
to develop new methodologies. Therefore, this study investigates the application of the
VFM formulated in 2D and 3D frameworks, respectively, 2D-VFM and 3D-VFM to the
identification of a hardening law up to large deformations. Even though only simulated
experiments are considered, the application of the VFM is investigated for a material with
different sample’s thickness to understand its influence.

The following section introduces the formulation of the VFM general to any 3D appli-
cation, as well as its numerical implementation. Then, the simulated tests and virtual
material used in this study are presented and analyzed. A section dedicated to the defi-
nition of virtual fields follows, and finally, the results of an analysis of predicted load and
identifications of a hardening law are presented.

**Virtual fields method**

Let us consider, in the Euclidean space, a body $\mathcal{B}$ subject to a general deformation process at a time instant $t$. According to the finite deformation theory [37], the position of its material points in the undeformed or reference configuration $\mathcal{B}_0$ (Lagrangian description) is given by $\mathbf{X}$, while in the deformed or current configuration $\mathcal{B}_t$ (Eulerian description)
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Description of the motion of a material body B in the Euclidean space defined by the global coordinate system ($\vec{x}, \vec{y}, \vec{z}$), its undeformed configuration, deformed configuration at a time instant t, and definition of a material point coordinate systems: ($\vec{i}, \vec{j}, \vec{k}$) is the global coordinate system defined in the undeformed configuration, ($\vec{1}, \vec{2}, \vec{3}$) is the local coordinate system defined in the deformed configuration, ($\vec{\Xi}, \vec{H}, \vec{\Pi}$) and ($\vec{\xi}, \vec{\eta}, \vec{\zeta}$) are the material coordinate systems, respectively, in the undeformed and deformed configurations. The angle $\theta$ represents the material orientation</figcaption></figure>

by **x**. The motion of each material point can be described as a function of **x** = $\chi(\mathbf{X}, t)$, mapping the position in the undeformed to the deformed configuration. The displacement field $\mathbf{u} = (u_x, u_y, u_z)$ which relates the position of each material point in the undeformed configuration and its position in the deformed configuration is given by

$$ \mathbf{u}(\mathbf{X}, t) = \chi(\mathbf{X}, t) - \mathbf{X} = \mathbf{x} - \mathbf{X}. \qquad (1) $$

Spatial derivatives can be used to calculate the deformation gradient F as

$$ \mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \frac{\partial \mathbf{u}}{\partial \mathbf{X}} + \mathbf{I}, \qquad (2) $$

where **I** is the second order identity tensor. The deformation gradient F provides a complete description of deformation, including stretch and rigid body rotations. The latter does not contribute to shape or size change, thus not contributing to the strain tensor. Therefore, it is necessary to separate the stretch from the rigid body rotations in the deformation gradient. Using the theorem of polar decomposition [38], F can be decomposed as

$$ \mathbf{F} = \mathbf{R}\mathbf{U} = \mathbf{V}\mathbf{R}, \qquad (3) $$

where **U** is the right stretch tensor, **V** is the left stretch tensor, and **R** is the orthogonal rotation tensor. A consequence of such mathematical description is that for every material point, a local coordinate system ($\vec{1}, \vec{2}, \vec{3}$), also known as corotational coordinate system, rotates during deformation, as schematically represented in Fig. 1.

The VFM is based on the principle of virtual work, which represents an integral form of mechanical equilibrium, written either in the undeformed or deformed configurations. Using the deformed configuration, for quasi-static conditions and neglecting body forces, the principle of virtual work can be given as

$$ \int_{\mathcal{B}_t} \sigma : \frac{\partial \mathbf{u}^*}{\partial \mathbf{x}} d\mathcal{B}_t + \int_{\partial \mathcal{B}_t} (\sigma \mathbf{n}) \cdot \mathbf{u}^* d\partial \mathcal{B}_t = 0, \qquad (4) $$
where $\sigma$ is the Cauchy stress tensor, $\mathbf{u}^*$ is the virtual displacements field, and $\mathbf{n}$ the outward vector of $\partial\mathcal{B}_t$. Alternatively, the principle of virtual work can be formulated in the undeformed configuration as

$$ \underbrace{\int_{\mathcal{B}_0} \mathbf{P} : \frac{\partial \mathbf{U}^*}{\partial \mathbf{X}} d\mathcal{B}_0}_{W_{\text{int}}} + \underbrace{\int_{\partial \mathcal{B}_0} (\mathbf{P}\mathbf{N}) \cdot \mathbf{U}^* d\partial \mathcal{B}_0}_{W_{\text{ext}}} = 0, \qquad (5) $$

where $\mathbf{P}$ is the first Piola-Kirchhoff stress tensor, the virtual displacements' field is now represented by $\mathbf{U}^*$, and $\mathbf{N}$ is the outward vector of $\partial\mathcal{B}_0$. The principle of virtual work expressed in the undeformed configuration is more practical for implementation since $\mathbf{U}^*$ (continuous and piecewise differentiable) remain unchanged and does not require updates to accommodate virtual boundary conditions [20,39].

The first Piola-Kirchhoff stress tensor $\mathbf{P}$ relates forces applied in the deformed configuration to volumes in the undeformed configuration, and is defined as

$$ \mathbf{P} = \det(\mathbf{F})\sigma\mathbf{F}^{-T}, \qquad (6) $$

where $\det(\mathbf{F})$ is the determinant of $\mathbf{F}$ and $\mathbf{F}^{-T}$ the transpose of its inverse. The Cauchy stress tensor $\sigma$ is calculated from the strain tensor $\varepsilon$ by means of a constitutive equation combined with a given set of material parameters $\xi$.

According to Eq. 5, the principle of virtual work states that, in equilibrium and for a continuous virtual field, the internal virtual work $W_{\text{int}}$ is equal to the external virtual work $W_{\text{ext}}$ in absolute value, which must be satisfied for any time instant. Let us now consider that the body $\mathcal{B}$ can be discretized by a given number of material points (similarly to a finite element mesh), then $W_{\text{int}}$ can be approximated by a discrete sum as

$$ W_{\text{int}}(\xi, \mathbf{U}^*, t) = \sum_{i=1}^{n_p} \left( \mathbf{P}_i(\xi, t) : \frac{\partial \mathbf{U}_i^*}{\partial \mathbf{X}_i} V_i \right), \qquad (7) $$

where $n_p$ represents the number of material points and $V_i$ is the volume covered by $i$-th material point. Due to the formulation of the principle of virtual work in the undeformed configuration, $V_i$ is constant for every time instant. Additionally, the external virtual work $W_{\text{ext}}$ can be simplified by selecting a virtual field that is constant along the boundary $\partial\mathcal{B}_0$ as

$$ \begin{aligned} W_{\text{ext}}(\mathbf{U}^*, t) &= \int_{\partial\mathcal{B}_0} (\mathbf{P}\mathbf{N}) \cdot \mathbf{U}^* d\partial\mathcal{B}_0 \\ &= \mathbf{U}^* \int_{\partial\mathcal{B}_0} (\mathbf{P}\mathbf{N}) d\partial\mathcal{B}_0 \\ &= \mathbf{U}^* \cdot \mathbf{L}(t), \end{aligned} \qquad (8) $$

where $\mathbf{L}$ is the resultant of the load acting on $\partial\mathcal{B}_t$, and recorded during the test. Because the load distribution on $\partial\mathcal{B}_0$ is usually unknown, this simplification is crucial for the VFM.

Finally, the VFM is used to identify the best set of material parameters by minimizing the difference between the internal and external virtual work, through an objective function
$\varphi$ written in a least-squares formulation as

$$ \varphi(\xi) = \sum_{t=1}^{n_t} \sum_{j=1}^{n_v} \left[ W_{\text{int}}(\xi, \mathbf{U}_j^*, t) - W_{\text{ext}}(\mathbf{U}_j^*, t) \right]^2, \qquad (9) $$

where $n_v$ is the number of virtual fields and $n_t$ the total number of time instants. Various time instants and different virtual fields can be used to enrich the objective function. Usually, an optimization algorithm is used to iteratively adjust the material parameters.

## Stress reconstruction from the strain fields

As previously established, the VFM requires, at each iteration, the calculation of $\mathbf{P}$, which is effectively dependent on the strain tensor. Therefore, it is convenient to use a measure of strain that is independent of rigid body rotations and dependent on stretch alone, given by

$$ \boldsymbol{\epsilon} = \ln \mathbf{V}, \qquad (10) $$

where $\mathbf{V}$ is defined in the global coordinate system ($\vec{i}, \vec{j}, \vec{k}$) in the undeformed configuration. However, constitutive equations are generally defined in the material coordinate system ($\vec{\xi}, \vec{\eta}, \vec{\zeta}$), aligned with the material orientation. Thus, the representation of the strain tensor in the global coordinate system $\boldsymbol{\epsilon}(\vec{i}, \vec{j}, \vec{k})$, must be expressed in the material coordinate system ($\vec{\xi}, \vec{\eta}, \vec{\zeta}$) according to

$$ \boldsymbol{\epsilon}(\vec{\xi}, \vec{\eta}, \vec{\zeta}) = \mathbf{R}_{\text{mat}}^{\text{T}} (\mathbf{R}^{\text{T}} \boldsymbol{\epsilon}(\vec{i}, \vec{j}, \vec{k}) \mathbf{R}) \mathbf{R}_{\text{mat}}, \qquad (11) $$

where $\mathbf{R}_{\text{mat}}$ is the material rotation tensor given by

$$ \mathbf{R}_{\text{mat}} = \begin{bmatrix} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{bmatrix}, \qquad (12) $$

which is responsible for rotating the local coordinate system to the material one as a function of the material orientation $\theta$, constant throughout deformation.

The projected strain tensor $\boldsymbol{\epsilon}(\vec{\xi}, \vec{\eta}, \vec{\zeta})$ is used to reconstruct the stress tensor in the material coordinate system $\boldsymbol{\sigma}(\vec{\xi}, \vec{\eta}, \vec{\zeta})$ via an implicit integration algorithm, similarly to the ones found in FEA software. However, many time instants should be used to guarantee the stability and convergence of the algorithm. Alternatively, other algorithms could be used, such as the direct method proposed by Rossi et al. [40], that provide a computationally faster reconstruction but cannot be incorporated in a FEA software.

Afterwards, $\boldsymbol{\sigma}(\vec{\xi}, \vec{\eta}, \vec{\zeta})$ should be expressed in the global coordinate system ($\vec{i}, \vec{j}, \vec{k}$) according to

$$ \boldsymbol{\sigma}(\vec{i}, \vec{j}, \vec{k}) = \mathbf{R}_{\text{mat}} (\mathbf{R} \boldsymbol{\sigma}(\vec{\xi}, \vec{\eta}, \vec{\zeta}) \mathbf{R}^{\text{T}}) \mathbf{R}_{\text{mat}}^{\text{T}}, \qquad (13) $$

and used in Eq. 5 to calculate the first Piola-Kirchhoff stress tensor $\mathbf{P}$ defined by Eq. 5.
The formulations presented up to here are general to any 3D application, thus correspond to the 3D-VFM. However, the VFM is commonly used with full-field measurements limited to the surface of the sample, corresponding to the 2D-VFM. These measurements only provide in-plane kinematic quantities, and often plane stress conditions are assumed, provided the sample's thickness is small relative to its width and height. Thus, the application is considered 2D, and out-of-plane stress components are negligible ($\sigma_{zz} = \sigma_{xz} = \sigma_{yz} = 0$), but the deformation gradient $\mathbf{F}$ is in reality fully 3D. It results that the calculation of $\mathbf{P}$ requires additional assumptions [39]. Assuming a thin specimen, the out-of-plane normal and shear components can be neglected ($F_{xz} = F_{yz} = F_{zx} = 0$), and the determinant of $\mathbf{F}$ is redefined as

$$ \det(\mathbf{F}) = F_{zz}(F_{xx}F_{yy} - F_{xy}F_{yx}), \qquad (14) $$

where $F_{xx}$, $F_{yy}$, $F_{xy}$, and $F_{yx}$ are the in-plane components of $\mathbf{F}$, and $F_{zz}$ is the out-of-plane component. The in-plane components are obtained from the displacement field. However, $F_{zz}$ is unknown, but it can be approximated by the out-of-plane strain component for each time instant $\varepsilon_{zz}(t)$ as

$$ F_{zz}(t) = 1 + \varepsilon_{zz}(t) = 1 + \int_{0}^{t} \dot{\varepsilon}_{zz} dt. \qquad (15) $$

Finally, assuming isotropic elastic behavior using Hooke's law and the isochoric behavior of plasticity, $\dot{\varepsilon}_{zz}$ can be determined as

$$ \dot{\varepsilon}_{zz} = -\frac{\nu}{1-\nu}(\dot{\varepsilon}_{xx}^{e} + \dot{\varepsilon}_{yy}^{e}) - (\dot{\varepsilon}_{xx}^{p} + \dot{\varepsilon}_{yy}^{p}), \qquad (16) $$

where $\dot{\varepsilon}_{xx}^{e}$ and $\dot{\varepsilon}_{yy}^{e}$ represent the components of elastic strain rate, and $\dot{\varepsilon}_{xx}^{p}$ and $\dot{\varepsilon}_{yy}^{p}$ the components of plastic strain rate.

## Numerical implementation

The VFM is numerically implemented in a code fully developed in Python programming language. In terms of computational efficiency, the code is developed with array-based calculations, instead of using loops over the material points. Using arrays to directly perform the operations, such as calculating the logarithmic strain of all material points, largely speeds up the computation and reduces differences into an implementation in Fortran. The small loss in computational speed is compensated in simplicity. Moreover, the array-based calculation easily allows the program to perform the operations in either a 2D or 3D framework.

As schematically presented in Fig. 2, the VFM requires the input of the geometry (coordinates of undeformed sconfiguration), displacement fields, evolution of load, and constitutive equations. It should be noted that this implementation uses a FE mesh and shape functions similar to a FEA software (such as Abaqus/Standard). In the pre-processing stage, the deformation gradient is calculated, and the logarithmic strain tensor is obtained from the polar decomposition of the deformation gradient. The latter operation is performed using a function available in the NumPy library for singular value decomposition [41]. The logarithmic strain tensor calculated in global coordinate system is rotated to the
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Flowchart of the VFM implementation</figcaption></figure>

material coordinate system in the undeformed configuration. Then, the virtual fields are generated manually.

The identification procedure uses the gradient-free Nelder-Mead (NM) algorithm included in the SciPy library [42,43]. This algorithm was chosen, in part, because it allows more control over the optimization procedure compared to the implementation of the Levenberg-Marquardt algorithm. For example, using the NM algorithm it is possible to call a function at the end of each iteration, which is not available for all algorithms of the SciPy library. This functionality is important to allow updates during the optimization procedure. Inside the identification procedure, the first and main operation is the calculation of the Cauchy stress tensor. Specifically, the Unified Material Modeling Driver for Plasticity (UMMDp) is integrated in the VFM to reconstruct the stress field through the backward-Euler scheme combined with an elastic predictor and plastic corrector [44]. The major advantages of using the UMMDp is to guarantee that the same model is used both in the identification procedure and future FE simulations, and to easily identify different constitutive models. A driver code for the UMMDp is developed that allows the communication with the main program through the F2PY framework [45]. This operation requires the most computational effort, because using an array-based calculation is impossible. The stress tensor is history dependent, so the calculation is performed iteratively for each material point and time instant. Afterwards, the Cauchy stress tensor is rotated to the global coordinate system, and the first Piola-Kirchhoff stress tensor is calculated. The internal and external virtual works are calculated, and finally, the objective function. Finally, the stopping criteria are evaluated, and the program either ends or updates the material parameters and continues the optimization procedure.
# Simulated tensile tests

## Specimen geometry and finite element model

As the aim is to give insights when 3D-VFM is necessary beyond necking in tension, this study considers, as a first step, only numerical simulations of a uniaxial tensile test. A specimen with varying cross-section is selected [33,46] and its geometry is presented in Fig. 3. Relatively to a classic uniaxial tensile specimen, this geometry has the advantage of generating a heterogeneous strain field over a large area, providing richer mechanical information for the identification of material parameters. Moreover, the varying cross-section will trigger necking in the center region of the specimen, facilitating the setup of acquisition cameras during real experiments.

The test is simulated with Abaqus/Standard [47], considering a fully 3D model of the test, without assuming symmetry conditions. The FE mesh is built using the 8-node brick elements with reduced integration (C3D8R), and considering three layers through the thickness. Ideally, more layers could have been used to describe more accurately the deformations along the thickness. However, it was found that increasing the number of layers also largely increased the computational time of numerical simulations and post-processing procedures. The FE mesh was generated with an in-plane global element size of 0.2 mm, resulting in a total of 197 955 elements. Nevertheless, the center region of the specimen is where deformation localizes, thus a volume of interest (VOI), represented by the colored region in Fig. 3, is defined and constituted by 97 307 elements. The in-plane mesh size was defined after performing a convergence study evaluating the maximum equivalent plastic strain, by considering only three layers of elements through the thickness. The selected mesh size provides the best compromise between convergence of the equivalent plastic strain values and the computational effort.

A Cartesian coordinate system is defined with the *x* axis along $d_w$, the *y* axis along $d_h$, and the *z* axis along $d_t$. The constants $d_w$ and $d_h$ represent, respectively, half of the width and length of the VOI, while $d_t$ is the sample's thickness. The origin of the coordinate system is located in the front surface of the specimen ($x = y = z = 0$), while $z = d_t$ defines the back surface. The numerical simulations are displacement-driven, with an imposed constant vertical displacement $u_y$ of 5 mm applied across the top edge, which is additionally fixed along the *x* and *z* axes to simulate the effect of the grip of a testing machine. The bottom edge of the model is fixed. The numerical simulations are performed at a constant time step, resulting in 200 time instants. Finally, a large strain formulation is adopted, and the constitutive model is implemented using the UMMDp.

## Virtual material

The material considered in this study is based on the real characteristics of a DP780 dual-phase steel with a nominal sample's thickness of 1.5 mm. This material was experimentally tested using quasi-static classical uniaxial tensile tests at different material orientations. The recorded longitudinal strain just before rupture was approximately 0.65. However, it should be emphasized that the goal of this study is not to characterize this material, but only to use it as a basis for the virtual material considered throughout. Multiple sample's thickness values of the same material are considered in this study, using the same description of the elastoplastic behavior.
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Geometry and finite element model of the specimen used in the numerical simulations of the simulated tensile test (dimensions in mm)</figcaption></figure>

<table><caption>**Table 1** Elastic properties of the material and reference parameters used in the Swift isotropic hardening law</caption><thead><tr><th colspan="2">Elasticity</th><th colspan="3">Swift hardening law</th></tr><tr><th><i>E</i> [GPa]</th><th><i>ν</i> [−]</th><th><i>K</i> [MPa]</th><th><i>σ</i><sub>0</sub> [MPa]</th><th><i>n</i> [−]</th></tr></thead><tbody><tr><td>200.0</td><td>0.3</td><td>1250.0</td><td>350.0</td><td>0.150</td></tr></tbody></table>

The elastic isotropic behavior of the material is described by Hooke's law. In plasticity, the yield condition $f(\boldsymbol{\sigma}, \bar{\varepsilon}^{\mathrm{P}})$, assuming isotropic hardening, can be expressed as

$$f(\boldsymbol{\sigma}, \bar{\varepsilon}^{\mathrm{P}}) = \bar{\sigma}(\boldsymbol{\sigma}) - \sigma_{\mathrm{y}}(\bar{\varepsilon}^{\mathrm{P}}) = 0, \qquad (17)$$

where $\bar{\sigma}(\boldsymbol{\sigma})$ is the equivalent stress and $\sigma_{\mathrm{y}}(\bar{\varepsilon}^{\mathrm{P}})$ is the yield stress as a function of the equivalent plastic strain $\bar{\varepsilon}^{\mathrm{P}}$. The elastoplastic behavior is described by the von Mises isotropic yield criterion and the Swift isotropic hardening law written as

$$\sigma_{\mathrm{y}}(\bar{\varepsilon}^{\mathrm{P}}) = K(\varepsilon_{0} + \bar{\varepsilon}^{\mathrm{P}})^{n}, \quad \text{with} \qquad (18)$$

$$\varepsilon_{0} = \left(\frac{\sigma_{0}}{K}\right)^{1/n}, \qquad (19)$$

where $K$ is the strength coefficient, $\varepsilon_0$ is the initial yield strain, $n$ is the strain-hardening exponent, and $\sigma_0$ is the initial yield stress. The latter model was selected for its simplicity and wide use in the scientific community. The elastic properties of the material and the Swift's law material parameters are provided in Table 1. The parameters of Swift law were obtained by fitting the experimental stress–strain curve of DP780. Additionally, in Fig. 4, the hardening behavior according to the Swift's law is presented up to 0.75 of equivalent plastic strain.

The material is defined with sample's thickness of 1.5, 2, 2.5, 3, 4, 6, 8, and 10 mm. The vertical displacement of 5 mm is used in all numerical simulations to have a direct comparison between samples. As a result, only the numerical model of the 1.5 mm thick sample is, approximately, modeled to reach rupture around the same magnitude of values as the real material. The equivalent plastic strain map for this sample, on the surface ($z = 0$) and on the middle plane ($z = d_t/2$), is shown in Fig. 5a at the last time instant corresponding to $u_y = 5$ mm. A large strain concentration is observed in the central region, indicating a severe necking just before rupture. For a smaller displacement of $u_y = 4$ mm (see Fig. 5a), it can be observed that the material was already in diffuse necking, though the values of equivalent plastic strain are much lower than at the last
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Virtual representation of the material strain-hardening behavior using the Swift isotropic hardening law and reference parameters</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5** Maps of equivalent plastic strain at time instants corresponding to vertical displacements of 2, 4, and 5 mm for samples' thickness of (a) 1.5 mm and (b) 10 mm. Maps on the surface ($z = 0$) and middle plane ($z = d_t/2$) are both represented</figcaption></figure>

time instant. This difference shows that large deformations are developed in the later stages of the test. On the other hand, the tests with large sample's thickness (such as 8 mm or 10 mm) do not even reach localized necking. It is observed that the necking is delayed as a function of the sample's thickness. This is observed in the equivalent plastic strain maps for the 10 mm thick sample presented in Fig. 5b. The values of equivalent plastic strain at the last time instant are only comparable to the ones observed at $u_y = 4$ mm for the 1.5 mm thick sample. The evolution of load as a function of the vertical displacement is also shown in Fig. 6 for all samples. As expected, the magnitude of load increases with the sample's thickness.

## Extraction of the displacement fields

The simulated tests are modeled in 3D, but the 2D-VFM uses information from the surface of the specimen, thus only in-plane kinematic fields are considered. For the 3D-VFM, the displacement field is extracted from the nodes of the FE mesh in the VOI of the simulated
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Evolution of the load as a function of the vertical displacement for the simulated tests with samples' thickness of (a) 1.5, 2, 2.5, and 3 mm, and (b) 4, 6, 8, and 10 mm</figcaption></figure>

tests. This mesh is already regularized, thus can be directly used as input to the 3D-VFM. For the 2D-VFM, the in-plane displacement field is extracted from the nodes on the surface of the 3D test, and a 2D FE mesh is generated by rearranging the elements' connectivity.

## Definition of the virtual fields

The main area of development and discussion around the VFM is probably the selection of virtual fields. Because its definition assigns weights to material points, it directly affects the objective function and can lead to different identification results. Nevertheless, one thing all methodologies agree is that the virtual fields must be kinematically admissible, continuous and differentiable over the entire domain [23].

In the framework of plasticity, two strategies can be distinguished according to the way virtual fields are generated. The most common strategy is based on manually defined virtual fields. This strategy usually employs polynomial or harmonic functions defined by the user and adapted to the geometry and boundary conditions of the application [31, 48, 49]. For this reason, this strategy highly depends on the user experience and understanding of the application, even though it is simpler from an implementation point of view. The other strategy is based on the automatic generation of virtual fields, and was developed to overcome the drawbacks of the first strategy. A major advantage is that it reduces the user dependence of the first strategy, making the method more suited for integration in parameter identification software. A first implementation, named optimized virtual fields, relied on the minimization of noise effects in the material parameters [50]. Later, Marek et al. [51] proposed the sensitivity-based virtual fields. This implementation relies on the sensitivity of the stress field to each material parameter, and it has been shown to outperform the manual virtual fields, particularly when dealing with experimental data [20, 28]. Nevertheless, it has been shown that increasing the number of manual virtual fields can lead to interesting results [52].

Considering that 3D-VFM is still less explored than 2D-VFM, only the manually defined virtual fields are employed in this study because of its simplicity and ease of implementation. The generic form of the virtual displacement fields U, for the general case of 3D, is given in the global coordinate system according to
<table><caption>**Table 2** Components of the virtual displacements and gradient of the virtual displacements for the three virtual fields</caption><thead><tr><th rowspan="2">Component</th><th colspan="3">Virtual field</th></tr><tr><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td colspan="4">Virtual displacements</td></tr><tr><td>U<sup>*</sup><sub>x</sub></td><td>0</td><td>1/π sin(π/x) cos(π/2d<sub>h</sub>) y</td><td>1/π sin(π/x) cos(π/2d<sub>h</sub>) y cos(π/2d<sub>t</sub>) z</td></tr><tr><td>U<sup>*</sup><sub>y</sub></td><td>y/d<sub>h</sub></td><td>1/π sin(π/x) cos(π/2d<sub>h</sub>) y</td><td>1/π sin(π/x) cos(π/2d<sub>h</sub>) y cos(π/2d<sub>t</sub>) z</td></tr><tr><td>U<sup>*</sup><sub>z</sub></td><td>0</td><td>0</td><td>1/π sin(π/x) cos(π/2d<sub>h</sub>) y sin(π/2d<sub>t</sub>) z</td></tr><tr><td colspan="4">Gradient of virtual displacements</td></tr><tr><td>∂U<sup>*</sup><sub>x</sub>/∂X</td><td>0</td><td>1/d<sub>w</sub> cos(π/x) cos(π/2d<sub>h</sub>) y</td><td>1/d<sub>w</sub> cos(π/x) cos(π/2d<sub>h</sub>) y cos(π/2d<sub>t</sub>) z</td></tr><tr><td>∂U<sup>*</sup><sub>x</sub>/∂Y</td><td>0</td><td>-1/2d<sub>h</sub> sin(π/x) sin(π/2d<sub>h</sub>) y</td><td>-1/2d<sub>h</sub> sin(π/x) sin(π/2d<sub>h</sub>) y cos(π/2d<sub>t</sub>) z</td></tr><tr><td>∂U<sup>*</sup><sub>y</sub>/∂Z</td><td>0</td><td>0</td><td>-1/2d<sub>t</sub> sin(π/x) cos(π/2d<sub>h</sub>) y sin(π/2d<sub>t</sub>) z</td></tr><tr><td>∂U<sup>*</sup><sub>y</sub>/∂X</td><td>0</td><td>1/d<sub>w</sub> cos(π/x) cos(π/2d<sub>h</sub>) y</td><td>1/d<sub>w</sub> cos(π/x) cos(π/2d<sub>h</sub>) y cos(π/2d<sub>t</sub>) z</td></tr><tr><td>∂U<sup>*</sup><sub>y</sub>/∂Y</td><td>1/d<sub>h</sub></td><td>-1/2d<sub>h</sub> sin(π/x) sin(π/2d<sub>h</sub>) y</td><td>-1/2d<sub>h</sub> sin(π/x) sin(π/2d<sub>h</sub>) y cos(π/2d<sub>t</sub>) z</td></tr><tr><td>∂U<sup>*</sup><sub>z</sub>/∂Z</td><td>0</td><td>0</td><td>-1/2d<sub>t</sub> sin(π/x) cos(π/2d<sub>h</sub>) y sin(π/2d<sub>t</sub>) z</td></tr><tr><td>∂U<sup>*</sup><sub>z</sub>/∂X</td><td>0</td><td>0</td><td>1/d<sub>w</sub> cos(π/x) cos(π/2d<sub>h</sub>) y sin(π/2d<sub>t</sub>) z</td></tr><tr><td>∂U<sup>*</sup><sub>z</sub>/∂Y</td><td>0</td><td>0</td><td>-1/2d<sub>h</sub> sin(π/x) sin(π/2d<sub>h</sub>) y sin(π/2d<sub>t</sub>) z</td></tr><tr><td>∂U<sup>*</sup><sub>z</sub>/∂Z</td><td>0</td><td>0</td><td>1/2d<sub>t</sub> sin(π/x) cos(π/2d<sub>h</sub>) y cos(π/2d<sub>t</sub>) z</td></tr></tbody></table>

$$ \mathbf{U}^* = \begin{bmatrix} U_x^* \\ U_y^* \\ U_z^* \end{bmatrix}, \qquad (20) $$

where $U_z^* = 0$ in the 2D-VFM, and the gradient of virtual displacement fields $\frac{\partial \mathbf{U}}{\partial X}$ by

$$ \frac{\partial \mathbf{U}^*}{\partial \mathbf{X}} = \begin{bmatrix} \frac{\partial U_x^*}{\partial X} & \frac{\partial U_x^*}{\partial Y} & \frac{\partial U_x^*}{\partial Z} \\ \frac{\partial U_y^*}{\partial X} & \frac{\partial U_y^*}{\partial Y} & \frac{\partial U_y^*}{\partial Z} \\ \frac{\partial U_z^*}{\partial X} & \frac{\partial U_z^*}{\partial Y} & \frac{\partial U_z^*}{\partial Z} \end{bmatrix}, \qquad (21) $$

where $\partial U_x^*/\partial Z = \partial U_y^*/\partial X = \partial U_z^*/\partial Z = \partial U_x^*/\partial Y = \partial U_z^*/\partial Z = 0$ in the 2D-VFM.

Three manual virtual fields are defined in this study, and their components are presented in Table 2. The first virtual field $\mathbf{U}_1^*$ is commonly used in various applications of the VFM because it includes the load's contribution [20,28,49]. The external virtual work is proportional to the load and is balanced by the internal virtual work calculated from the stress component $\sigma_{yy}$. In theory, it is the only virtual field that gives a non-zero value in the calculation of the internal virtual work. Moreover, the resulting gradient presents a constant spatial distribution across all VOI.

The second and third virtual fields, respectively, denoted by $\mathbf{U}_2^*$ and $\mathbf{U}_3^*$, present non-constant spatial distributions, with the exception of the z component in the second virtual field. The varying spatial distributions lead to different weights assigned to the stress components in the calculation of the internal virtual work. These are null on the boundary
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Maps of virtual displacement fields represented for three virtual fields (a) $\mathbf{U}_1^*$, (b) $\mathbf{U}_2^*$, and (c) $\mathbf{U}_3^*$. The maps of $\mathbf{U}_1^*$ are represented on the front ($z = 0$) and back ($z = d_t$) surfaces, while the maps of $\mathbf{U}_1^*$ and $\mathbf{U}_3^*$ are representative of the whole VOI</figcaption></figure>

of the applied load, which gives a zero value in the calculation of external virtual work. The second virtual field was adopted from other studies using 2D-VFM, where its use has proven to be successful [20, 31, 39]. It only considers the in-plane components of the stress tensor, neglecting the out-of-plane components. Also, $\mathbf{U}_2^*$ does not depend on $z$, leading to a constant spatial distribution across the thickness. It should be noted that the extension of this virtual field to a 3D formulation is, to the best of the author's knowledge, a novelty.

Lastly, $\mathbf{U}_3^*$ assigns weight to all the stress tensor components and also generates a non-constant spatial distribution across the thickness. The third virtual field can be considered an extension of the second to 3D-VFM, because it simplifies to the latter when $z = 0$. However, it should be noted that this virtual field would not be useful if the origin of the coordinate system was positioned in the middle of the specimen, because the kinematic fields are symmetrical across the thickness. Thus, leading to a combination of symmetrical stress fields with symmetrical virtual fields across the thickness, resulting in a null value of the internal virtual work. The former is only true because this study only considers simulated tests, while in real experiments unsymmetrical kinematic fields across the thickness could be expected. A graphical representation of the three virtual fields displacement maps is shown in Fig. 7. The presented maps of $\mathbf{U}_1^*$ and $\mathbf{U}_2^*$ are representative of the whole VOI, while for $\mathbf{U}_3^*$ are presented the maps at the front ($z = 0$) and back ($z = d_t$) surfaces.
# Results

## Analysis of the predicted load

Because this study only considers numerical data, the reference parameters and reference test data are known. By taking advantage of that, it is possible to analyze the predicted load using the 2D-VFM and the 3D-VFM, and compare with the reference load obtained from the numerical simulations. In particular, the evolution of load can be extracted from numerical simulations and used in the external virtual work term, as shown in Eq. 8, provided a constant virtual field over the boundary is used. Also, considering the test is uniaxially loaded along the y axis, the evolution of load can be simplified to its y component, denoted by $L^{\text{ref}}(t)$. Using the first virtual field $U_1^*$, the external virtual work $W_{\text{ext}}(U_1^*, t)$ can be simplified to

$$ W_{\text{ext}}(U_1^*, t) = \frac{y}{d_h} \cdot L^{\text{ref}}(t) \qquad (22) \\ = L^{\text{ref}}(t), $$

where $U_1^* = y/d_h = 1$ on the boundary of applied load. Consequently, the evolution of predicted load $L(\xi, t)$ can be calculated from the internal virtual work $W_{\text{int}}(\xi, U_1^*, t)$ as

$$ \begin{aligned} L(\xi, t) &= \frac{1}{2} W_{\text{int}}(\xi, U_1^*, t) \\ &= \frac{1}{2} \int_{B_0} \mathbf{P}(\xi, t) : \frac{\partial U_1^*}{\partial \mathbf{X}} d\mathcal{B}_0 \qquad (23) \\ &= \frac{1}{2} \int_{B_0} P_{yy}(\xi, t) \cdot \frac{1}{d_h} d\mathcal{B}_0, \end{aligned} $$

The constant $1/2$ is necessary to account for the reaction load on the fixed boundary.

The relative error $\delta_{\text{re}}(\xi, t)$ between predicted and reference load at a given time instant $t$ given by

$$ \delta_{\text{re}}(\xi, t) = 100 \times \frac{|L(\xi, t) - L^{\text{ref}}(t)|}{L^{\text{ref}}(t)}, \qquad (24) $$

is used as the indicator to compare the results using the 2D-VFM and the 3D-VFM.

In Fig. 8, the evolution of the relative error is presented for all samples. It is observed that using the 3D-VFM, the maximum relative error is lower than 0.03% for all samples, and tends to decrease with the increase in sample's thickness. Thus, it can be considered that the relative error of 3D-VFM is negligible, and that the evolution of load is successfully predicted. However, the same is not observed using the 2D-VFM, where the relative error is not negligible. Overall, the relative error increases with the increase in sample's thickness. For 1.5 mm, the evolution of relative error is close to zero up to 4 mm of displacement, corresponding to the maximum load, where the error rapidly increases up to 1%. It was shown before that the time instant corresponding to $u_y = 4$ mm is when the equivalent plastic strain values rapidly increase. Increasing the sample's thickness, an increase in the relative error is observed sooner in the test, even though the same displacement corresponds to lower deformations in the thicker tests. For the thicker sample of 10 mm, the relative error is as high as 6% at the last time instant, though it can be expected that it would be higher for later stages of the test leading to rupture. In general, the increase in relative error is related to the sample's thickness.
<figure><img src="image_11.png" /><figcaption>**Fig. 8** Evolution of the relative error for predicted load using the 2D-VFM and the 3D-VFM with the reference material parameters, for tests with thickness of (a) 1.5, 2, 2.5, and 3 mm, and (b) 4, 6, 8, and 10 mm</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9** Location of the material points $P_S$ and $P_M$, represented on the finite element mesh of the test with 1.5 mm thick sample</figcaption></figure>

To better understand the origin of these deviations, the evolution of some stress components is investigated for two material points, one in the center of the front surface, denoted $P_S$, and another in the center of the middle plane, denoted $P_M$, from the FE results. The exact location of these two material points in the finite element mesh of all tests is shown in Fig. 9, for the 1.5 mm thick sample as an example.

This analysis is only performed for the samples' thickness of 1.5, 3, and 10 mm, as presented in Fig. 10. It can be seen that there is a dominant stress component $\sigma_{yy}$ associated with the loading direction, which is expected from a uniaxial tensile test. More interestingly, it can be seen for the 1.5 mm thick sample that the time instant where the error
<figure><img src="image_13.png" /><figcaption>Fig. 10 Evolution of the stress components $\sigma_{xx}$, $\sigma_{yy}$, and $\sigma_{zz}$ for the material points $P_S$ and $P_M$, respectively, located on the center of the surface and in the center of the middle plane, for samples' thickness of (a) 1.5 mm, (b) 3 mm, and (c) 10 mm</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11 Evolution of ratio of stress components $\sigma_{xx}$, $\sigma_{yy}$, and $\sigma_{zz}$ between material points $P_S$ and $P_M$, for samples' thickness of (a) 1.5 mm, (b) 3 mm, and (c) 10 mm</figcaption></figure>

rapidly increases is consistent with the appearance of a non-zero stress component $\sigma_{zz}$ (see Fig. 10a). In fact, this is also consistent with the onset of localized necking, as previously observed in Fig. 4. For larger sample's thickness, the development of $\sigma_{zz}$ is not as sudden and intense because the localized necking has not been reached, but more as a constant increase from the initial time instant, in particular for the 10 mm thick sample (see Fig. 10c). Moreover, for the three samples, there is a difference in the stress values between the two material points.

Additionally, an analysis of the evolution of the ratio of stress components is also performed between the two material points. Its formulation is given by

$$R_{ii} = \frac{\sigma_{ii}(P_M)}{\sigma_{ii}(P_S)} \quad \text{with } i = x, y, z. \qquad (25)$$

This analysis shows that the magnitude of $\sigma_{zz}$ in $P_M$ is two times higher than in $P_S$ (see Fig. 11). The same behavior is observed for all samples, as shown in Table 3, where the maximum values of the ratio of stress are presented. It is also worth noting the tendency of $R_{xx}$ to increase for higher sample's thickness. Since the 2D-VFM only considers the in-plane stress components and the surface material points in the calculation of the internal virtual work, it is not able to capture these behaviors, leading to the errors in the calculation of the internal virtual work, and consequently, errors in the predicted load.

Although errors will exist when using the 2D-VFM, it is important to understand if it is an acceptable error. For example, a relative error of 1% (1.5 mm) can be considered acceptable for an application not requiring high accuracy. As such, it can also be considered that for sample's thickness up to 1.5 mm, using the 2D-VFM is acceptable. However, for
<table><caption>**Table 3** Maximum values of the ratio of stress components $\sigma_{xx}$, $\sigma_{yy}$, and $\sigma_{zz}$ between material points $P_S$ and $P_M$, for all samples</caption><thead><tr><th rowspan="2">Ratio of stress [−]</th><th colspan="8">Sample thickness [mm]</th></tr><tr><th>1.5</th><th>2</th><th>2.5</th><th>3</th><th>4</th><th>6</th><th>8</th><th>10</th></tr></thead><tbody><tr><th>R<sub>xx</sub></th><td>1.15</td><td>1.16</td><td>1.18</td><td>1.20</td><td>1.25</td><td>1.37</td><td>1.54</td><td>1.73</td></tr><tr><th>R<sub>yy</sub></th><td>1.07</td><td>1.04</td><td>1.04</td><td>1.03</td><td>1.03</td><td>1.04</td><td>1.05</td><td>1.06</td></tr><tr><th>R<sub>zz</sub></th><td>2.14</td><td>2.04</td><td>2.02</td><td>2.02</td><td>2.01</td><td>2.01</td><td>2.00</td><td>2.00</td></tr></tbody></table>

<table><caption>**Table 4** Results of the identification starting from the reference material parameters and using the 2D-VFM and the 3D-VFM, for tests with sample's thickness of 1.5, 3, and 10 mm. Only the first virtual field is employed in the identification</caption><thead><tr><th rowspan="2">Sample thickness</th><th colspan="3">Parameter</th><th colspan="2">Objective function</th></tr><tr><th>K [MPa]</th><th>σ<sub>0</sub> [MPa]</th><th>n [−]</th><th>Initial [−]</th><th>Final [−]</th></tr></thead><tbody><tr><td>initial</td><td>1250.000</td><td>350.000</td><td>0.150</td><td>–</td><td>–</td></tr><tr><td>2D-VFM</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1.5 mm</td><td>1240.609</td><td>340.502</td><td>0.148</td><td>3.00d6</td><td>2.10d6</td></tr><tr><td>3 mm</td><td>1217.838</td><td>318.345</td><td>0.143</td><td>9.10d7</td><td>3.30d7</td></tr><tr><td>10 mm</td><td>1137.783</td><td>5.049</td><td>0.128</td><td>1.58d10</td><td>2.40d9</td></tr><tr><td>3D-VFM</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1.5 mm</td><td>1250.591</td><td>350.365</td><td>0.150</td><td>3.31d3</td><td>8.20d2</td></tr><tr><td>3 mm</td><td>1250.463</td><td>350.255</td><td>0.150</td><td>7.76d3</td><td>8.30d2</td></tr><tr><td>10 mm</td><td>1249.936</td><td>349.863</td><td>0.150</td><td>7.11d3</td><td>5.40d3</td></tr></tbody></table>

larger sample's thickness, a larger deviation in the internal virtual work should be expected, leading to larger errors in the predicted load.

## Identification of the Swift hardening law

The developed VFM code in Python programming language is applied to identify the material parameters of the Swift's law, with $\xi = (K, \sigma_0, n)$. The adopted termination criteria of the optimization procedure are set to a maximum number of iterations of 500, or when the absolute difference between parameters in consecutive iterations is lower than or equal to $1 \times 10^{-8}$. The identification is carried out for only the three samples' thickness of 1.5, 3, and 10 mm.

In a first part, the identification is carried out using both the 2D-VFM and the 3D-VFM. Only the first virtual field is used here. The three material parameters of Swift's law are identified, and the starting guess is equal to the reference parameters (see Table 1). This way, it will be easier to verify if the reference parameters are the global minimum solution used by numerical simulations. It is expected that the optimization is easily carried out, without taking many iterations or stopping at local minima.

The results of this identification are summarized in Table 4, where it is observed that using the 3D-VFM, the final material parameters are identical to the reference ones. In fact, considering the objective function is not normalized, the initial and final values are of low magnitude. These results reinforce the previous analysis, where very low errors were obtained for the predicted load using the 3D-VFM.

However, the results using the 2D-VFM are completely different. By comparison, the initial objective function values are higher than those observed using the 3D-VFM, which directly indicates a worse prediction of the load using the reference material parameters. For the 1.5 mm thick sample, the identified material parameters are not that far from the
<figure><img src="image_15.png" /><figcaption>**Fig. 12** Evolution of the relative error for predicted load using the reference and the identified material parameters using the 2D-VFM, for tests with samples' thickness of (a) 1.5 mm, (b) 3 mm, and (c) 10.0 mm</figcaption></figure>

reference ones, even though it directly represents that these are suitable parameters. In fact, the objective function values have not largely improved, and its values are still much higher than those observed using the 3D-VFM. Moreover, the 2D-VFM uses a third of the material points considered in the 3D-VFM, thus proportionally the objective function values using 2D-VFM would be even higher. For the 3, and 10 mm thick samples, the identified parameters are far from the reference parameters. In particular, the identified value of $\sigma_0$ for the 10 mm thick sample is not even physically acceptable. To avoid unfeasible solutions, bounds to the parameters can be imposed.

In Fig. 12, the evolution of the relative error for predicted load (see Eq. 24) is represented using both the reference and identified material parameters from Table 4 using the 2D-VFM. For the 1.5 mm thick sample, it is observed that the predicted load with the identified parameters has improved in the range of necking (from approximately 4 mm of displacement). In spite of that, the relative error increases in early stages of the test, where the error was close to zero using the reference parameters. The same phenomenon is observed for the 3, and 10 mm thick samples with even more emphasis. It can be concluded that the identification using the 2D-VFM results in identified material parameters that are not representative of the materials' hardening behavior when applied to samples with large thickness.

To analyze the influence of each virtual field (see Table 2), more identifications are performed, but only using the 3D-VFM. In the identifications, various combinations of virtual fields are used. Three identifications evaluate the influence of each virtual field individually, and three more are performed by combining the first with the second virtual field, the second with the third, and a final identification combining the three of them. Again, the identifications are performed only for the 1.5, 3, and 10 mm thick samples. However, the initial material parameters are different from the reference ones to investigate the effect on the efficiency of convergence. The selected initial material parameters correspond to 1.5 times the reference ones.

Overall, the identification results are good, as summarized in Table 5. A common result between all identifications is the exact identification of the parameter *n*, while for the parameters *K* and $\sigma_0$ the exact values are not identified, but almost identical. The values of the identified parameters are all within 1% of the reference ones, which is more than an acceptable result. Nevertheless, for the sake of comparing the influence of different virtual fields, it can be said that the best results are achieved using the first virtual field, either alone or in combination with others. In fact, this virtual field is the only one that leads to a non-
<table><caption>**Table 5** Results of the identification using the 3D-VFM and combining different virtual fields, for tests with samples' thickness of 1.5, 3, and 10 mm</caption><thead><tr><th rowspan="2">Sample thickness</th><th rowspan="2">Virtual fields</th><th colspan="3">Parameter</th><th rowspan="2">Iter.</th><th colspan="2">Objective function</th></tr><tr><th>K [MPa]</th><th>σ₀ [MPa]</th><th>n [−]</th><th>Initial [−]</th><th>Final [−]</th></tr></thead><tbody><tr><td rowspan="6">1.5 mm</td><td>initial</td><td>1875.000</td><td>525.000</td><td>0.225</td><td>–</td><td>–</td><td>–</td></tr><tr><td>1</td><td>1250.591</td><td>350.365</td><td>0.150</td><td>217</td><td>7.64e9</td><td>8.21e2</td></tr><tr><td>2</td><td>1247.688</td><td>348.374</td><td>0.150</td><td>226</td><td>1.97e8</td><td>9.68e2</td></tr><tr><td>3</td><td>1247.355</td><td>348.439</td><td>0.150</td><td>217</td><td>8.25e7</td><td>8.72e2</td></tr><tr><td>1 + 2</td><td>1250.454</td><td>350.240</td><td>0.150</td><td>218</td><td>7.89e9</td><td>3.71e3</td></tr><tr><td>1 + 3</td><td>1250.525</td><td>350.307</td><td>0.150</td><td>240</td><td>6.73e9</td><td>2.70e3</td></tr><tr><td rowspan="6">3 mm</td><td>1 + 2 + 3</td><td>1250.394</td><td>350.187</td><td>0.150</td><td>232</td><td>6.93e9</td><td>5.52e3</td></tr><tr><td>1</td><td>1250.463</td><td>350.255</td><td>0.150</td><td>234</td><td>3.06e10</td><td>8.32e2</td></tr><tr><td>2</td><td>1248.317</td><td>348.930</td><td>0.150</td><td>208</td><td>9.13e8</td><td>7.28e2</td></tr><tr><td>3</td><td>1248.670</td><td>349.627</td><td>0.150</td><td>217</td><td>3.80e8</td><td>3.88e1</td></tr><tr><td>1 + 2</td><td>1250.357</td><td>350.160</td><td>0.150</td><td>216</td><td>3.15e10</td><td>6.45e3</td></tr><tr><td>1 + 3</td><td>1250.423</td><td>350.223</td><td>0.150</td><td>221</td><td>3.10e10</td><td>2.56e3</td></tr><tr><td rowspan="6">10 mm</td><td>1 + 2 + 3</td><td>1250.321</td><td>350.132</td><td>0.150</td><td>220</td><td>3.19e10</td><td>8.04e3</td></tr><tr><td>1</td><td>1249.936</td><td>349.863</td><td>0.150</td><td>240</td><td>3.51e11</td><td>5.39e3</td></tr><tr><td>2</td><td>1246.999</td><td>347.853</td><td>0.150</td><td>210</td><td>9.40e09</td><td>1.67e4</td></tr><tr><td>3</td><td>1249.914</td><td>349.843</td><td>0.150</td><td>275</td><td>3.09e11</td><td>1.09e4</td></tr><tr><td>1 + 2</td><td>1249.792</td><td>349.726</td><td>0.150</td><td>248</td><td>3.62e11</td><td>1.18e5</td></tr><tr><td>1 + 3</td><td>1249.914</td><td>349.843</td><td>0.150</td><td>275</td><td>3.09e11</td><td>1.09e4</td></tr><tr><td></td><td>1 + 2 + 3</td><td>1249.773</td><td>349.710</td><td>0.150</td><td>260</td><td>3.19e11</td><td>1.22e5</td></tr></tbody></table>

zero value in external virtual work, including valuable information to the identification of hardening material parameters. When using either the second or third virtual fields alone, the identified parameters K and σ₀ tend to not reach the exact solution. Nevertheless, this does not imply that these virtual fields are unsuitable; rather, it reinforces the idea that manually defining virtual fields and increasing their number leads to more accurate and stable results [52]. Interestingly, for the 10 mm thick sample, the results using the second virtual field are slightly worst compared to the results of the 1.5, and 3 mm thick samples. This finding suggests that when out-of-plane stress components are not negligible, using virtual fields similar to the second leaves out important information that could be captured in the objective function. In terms of computational efficiency, the number of iterations required for the convergence of the identification procedure does not suggest any major differences between the combination of virtual fields.

## Conclusions

The purpose of this study is to better understand the limitations of using the VFM in a 2D framework compared to a 3D framework. The study only considers data from numerical simulations, and the methodology was applied to the identification of a hardening law, beyond necking and up to rupture. The simulated tests are performed with a uniaxial tensile test of varying cross-sections, generating a heterogeneous strain field over a large area. Moreover, various samples of different thickness are considered to understand its influence on the use of the VFM. Although these samples present different thickness, the material's hardening behavior is based on the real behavior of a dual-phase steel with 1.5 mm of thickness.
Three manual virtual fields are implemented in this study, because of its simplicity and lack of studies in a 3D framework. The first virtual field is commonly used in the VFM, because it is the only one that includes the external virtual work in the calculation of the objective function. The second virtual field only assigns weight to in-plane stress components. It is also widely used in the VFM in a 2D framework. The third, and last, virtual field is an extension of the second to a 3D framework, and can be considered a novel contribution of this study. The advantage of the latter is that it assigns non-zero weights to all stress components and is not constant along the thickness.

In a first part, the evolution of the predicted load was analyzed using the 2D-VFM and the 3D-VFM, for all samples. The relative error between predicted and reference load from numerical simulations was used as an indicator. The results showed that the relative error was negligible when using the 3D-VFM, successfully predicting the evolution of the load. On the contrary, using the 2D-VFM, the relative error was no longer negligible, ranging from 1% for the 1.5 mm thick sample up to 6% for the 10 mm thick sample, at the last time instant. With this analysis, it has been shown that the relative error increases proportionally to the increase of the sample's thickness. Moreover, the source of the error can be related to the development of non-negligible out-of-plane stress components not considered in the 2D-VFM, as well as considerable differences in the stress values between the surface and middle of the specimen.

In a second part, the identification of the Swift hardening law material parameters was investigated. Initially, the identification was carried out using both VFM frameworks, only one virtual field, and starting from the reference material parameters. The results have shown that the reference material parameters correspond to the global optimum solution using the 3D-VFM, but using the 2D-VFM led to solutions diverging from the reference, and one unfeasible solution was also found. A verification of the predicted load using the identified parameters has shown that these do not lead to better results, but rather improve the prediction at some stages at the expense of others. Lastly, the three virtual fields were employed in the identification of the material parameters using the 3D-VFM. This time, a starting solution far away from the reference was chosen to analyze the convergence efficiency using different combinations of virtual fields. Nevertheless, few differences were observed between identifications, with all reaching solutions within 1% of the reference values or in some cases identical.

Finally, it can be said that this study performed a pioneering analysis of the limitations of the VFM in a 2D framework compared to a 3D framework. However, this was only a preliminary numerical study that can be enlarged and extended to materials developing pronounced out-of-plane stress fields.

## Author contributions

Not applicable.

## Funding Information

This project has received funding from the Research Fund for Coal and Steel under Grant agreement No. 888153. The authors also gratefully acknowledge the financial support of the Portuguese Foundation for Science and Technology (FCT) under the projects UIDB/00481/2020 and UIDP/00481/2020, and CENTRO-01-0145-FEDER-022083 by UE/FEDER through the programs CENTRO 2020 and COMPETE 2020. The authors also would like to acknowledge the RÂgion Bretagne for its financial support and BPI France under project EXPRESSo. Miguel G. Oliveira is grateful to the FCT for the PhD Grant SFRH/BD/143665/2019.
**Availability of data materials** The datasets used and/or analysed during the current study are available from the
corresponding author on reasonable request.

**Code availability** Not applicable.

**Material availability** Not applicable.

**Declarations**

**Ethics approval and consent to participate**

Not applicable.

**Consent for publication**

Not applicable.

**Competing interests**

The authors declare that they have no Conflict of interest.

Received: 27 May 2024 Accepted: 20 March 2025

Published online: 28 April 2025

**References**

1. Hill R. On discontinuous plastic states, with special reference to localized necking in thin sheets. J Mech Phys Solids. 1952;1(1):19–30. https://doi.org/10.1016/0022-5096(52)90003-3.

2. Hutchinson JW, Miles JP. Bifurcation analysis of the onset of necking in an elastic/plastic cylinder under uniaxial tension. J Mech Phys Solids. 1974;22(1):61–71. https://doi.org/10.1016/0022-5096(74)90014-3.

3. Nichols FA. Plastic instabilities and uniaxial tensile ductilities. Acta Metall. 1980;28(6):663–73. https://doi.org/10.1016/0001-6160(80)90144-3.

4. Tvergaard V. Necking in tensile bars with rectangular cross-section. Comput Methods Appl Mech Eng. 1993;103:273–90. https://doi.org/10.1016/0045-7825(93)90049-4.

5. Tong W, Zhang N. An experimental investigation of necking in thin sheets. In: Proceedings of the ASME 2001 International Mechanical Engineering Congress and Exposition. 2001;111-16:231–238.

6. Dunand M, Mohr D. Hybrid experimental-numerical analysis of basic ductile fracture experiments for sheet metals. Int J Solids Struct. 2010;47(9):1130–43. https://doi.org/10.1016/j.ijsolstr.2009.12.011.

7. Tarfid N, Kyriakides S. Determination of anisotropy and material hardening for aluminum sheet metal. Int J Solids Struct. 2012;49(25):3496–506. https://doi.org/10.1016/j.ijsolstr.2012.01.011.

8. Tu T, Ren X, He J, Zhang Z. Stress-strain curves of metallic materials and post-necking strain hardening characterization: a review. Fatigue Fracture Eng Mater Struct. 2020;43(1):3–19.

9. Bridgman PW. Studies in large plastic flow and fracture: with special emphasis on the effects of hydrostatic pressure. McGraw Hill; 1952.

10. Zhang ZL, Hauge M, Ødegård J, Thaulow C. Determining material true stress-strain curve from tensile specimens with rectangular cross-section. Int J Solids Struct. 1999;36(23):3497–516. https://doi.org/10.1016/s0020-7683(98)00153-x.

11. Mirone G. A new model for the elastoplastic characterization and the stress-strain determination on the necking section of a tensile specimen. Int J Solids Struct. 2004;41(13):3545–64. https://doi.org/10.1016/j.ijsolstr.2004.02.011.

12. Koc P, Štok B. Computer-aided identification of the yield curve of a sheet metal after onset of necking. Comput Mater Sci. 2004;31(1):155–68. https://doi.org/10.1016/j.commatsci.2004.02.004.

13. Tao H, Zhang N, Tong W. An iterative procedure for determining effective stress-strain curves of sheet metals. Int J Mech Mater Des. 2009;5(1):13–27. https://doi.org/10.1007/s10999-008-9082-2.

14. Pottier T, Toussaint F, Vacher P. Contribution of heterogeneous strain field measurements and boundary conditions modelling in inverse identification of material parameters. Eur J Mech A Solids. 2011;30(3):373–82. https://doi.org/10.1016/j.euromechsol.2010.10.001.

15. Kajberg J, Lindkvist G. Characterisation of materials subjected to large strains by inverse modelling based on in-plane displacement fields. Int J Solids Struct. 2004;41(13):3439–59. https://doi.org/10.1016/j.ijsolstr.2004.02.021.

16. Denys K, Coppieters S, Seefeldt M, Debruyne D. Multi-DIC setup for the identification of a 3D anisotropic yield surface of thick high strength steel using a double perforated specimen. Mech Mater. 2016;100:96–108. https://doi.org/10.1016/j.mechmat.2016.06.011.

17. Kacem A, Laurent H, Thuillier S. Influence of experimental boundary conditions on the calibration of a ductile fracture criterion. Eng Fract Mech. 2021;248: 107686. https://doi.org/10.1016/j.engfracmech.2021.107686.

18. Pradeau A, Thuillier S, Yoon JW. Prediction of failure in bending of an aluminium sheet alloy. Int J Mech Sci. 2016;119:23–35. https://doi.org/10.1016/j.ijmecsci.2016.09.033.

19. Rossi M, Lattanzi A, Barlat F. A general linear method to evaluate the hardening behaviour of metals at large strain with full-field measurements. Strain. 2018;54(3):12265. https://doi.org/10.1111/str.12265.

20. Lattanzi A, Barlat F, Pierron F, Marek A, Rossi M. Inverse identification strategies for the characterization of transformation-based anisotropic plasticity models with the non-linear VFM. Int J Mech Sci. 2020;173: 105422. https://doi.org/10.1016/j.ijmecsci.2020.105422.

21. Coppieters S, Cooreman S, Sol H, Van Houtte P, Debruyne D. Identification of the post-necking hardening behaviour of sheet metal by comparison of the internal and external work in the necking zone. J Mater Process Technol. 2011;211(3):545–52. https://doi.org/10.1016/j.jmatprotec.2010.11.015.

22. Coppieters S, Kuwabara T. Identification of post-necking hardening phenomena in ductile sheet metal. Exp Mech. 2014;54(8):1355–71. https://doi.org/10.1007/s11340-014-9900-4.
23. Pierron F, Grédiac M. The virtual fields method: extracting constitutive mechanical parameters from full-field deformation measurements. Springer, 2012.

24. Avril S, Huntley JM, Pierron F, Steele DD. 3D Heterogeneous stiffness reconstruction using MRI and the virtual fields method. Exp Mech. 2008;48(4):479-94. https://doi.org/10.1007/s11340-008-9128-2.

25. Notta-Cuvier D, Langrand N, Lauro F, Markiewicz E. An innovative procedure for characterising a coupled elastoplastic damage model of behaviour using the virtual fields method. Int J Solids Struct. 2015;69-70:415-27. https://doi.org/10.1016/j.ijsolstr.2015.05.009.

26. Jones E, Karlson KN, Reu PL. Investigation of assumptions and approximations in the virtual fields method for a viscoplastic material model. Strain. 2019;55(4):e12309. https://doi.org/10.1111/str.12309.

27. Martins J, Andrade-Campos A, Thuillier S. Calibration of anisotropic plasticity models using a biaxial test and the virtual fields method. Int J Solids Struct. 2019;172-173:21-37. https://doi.org/10.1016/j.ijsolstr.2019.05.019.

28. Marek A, Davis FM, Kim JH, Pierron F. Experimental validation of the sensitivity-based virtual fields for identification of anisotropic plasticity models. Exp Mech. 2020;60(5):639-64. https://doi.org/10.1007/s11340-019-00575-3.

29. Fu J, Zhu K, Nie X, Tang Y, Yang Z, Qi L. Inertia-based identification of elastic anisotropic properties for materials undergoing dynamic loadings using the virtual fields method and heterogeneous impact tests. Mater Design. 2021;203: 109594. https://doi.org/10.1016/j.matdes.2021.109594.

30. Rossi M, Lattanzi A, Morichelli L. Integrated thermomechanical characterization of metals using the virtual fields method. Key Eng Mater. 2022;926:2222-7. https://doi.org/10.4028/p-t8ksaq.

31. Rossi M, Pierron F, Štamborská M. Application of the virtual fields method to large strain anisotropic plasticity. Int J Solids Struct. 2016;97-98:322-35. https://doi.org/10.1016/j.ijsolstr.2016.07.015.

32. Coppieters S, Kim JH, Denys K, Cooreman S, Debruyne D. On complete solutions for the problem of diffuse necking in sheet metal. Procedia Eng. 2017;207:2012-7. https://doi.org/10.1016/j.proeng.2017.10.1060.

33. Kim JH, Serpantié A, Barlat F, Pierron F, Lee MG. Characterization of the post-necking strain hardening behavior using the virtual fields method. Int J Solids Struct. 2013;50(24):3829-42. https://doi.org/10.1016/j.ijsolstr.2013.07.018.

34. Rossi M, Cortese L, Genovese K, Lattanzi A, Nalli F, Pierron F. Evaluation of volume deformation from surface DIC measurement. Exp Mech. 2018;58(7):1181-94. https://doi.org/10.1007/s11340-018-0409-0.

35. Oliveira MG, Lattanzi A, Thuillier S, Andrade-Campos A, Rossi M. Error estimation of the 3D reconstruction of a displacement field from DIC measurements. Key Eng Mater. 2022;926:2129-36. https://doi.org/10.4028/p-abe455.

36. Rossi M, Pierron F. Identification of plastic constitutive parameters at large deformations from three dimensional displacement fields. Comput Mech. 2012;49(1):53-71. https://doi.org/10.1007/s00466-011-0627-0.

37. Dunne F, Petrinic N. Introduction to Computational Plasticity. Oxford University Press; 2005.

38. de Souza Neto EA, Perić D, Owen D. Computational methods for plasticity: theory and applications. John Wiley & Sons; 2008.

39. Marek A, Davis FM, Rossi M, Pierron F. Extension of the sensitivity-based virtual fields to large deformation anisotropic plasticity. Int J Mater Form. 2019;12(3):457-76. https://doi.org/10.1007/s12289-018-1428-1.

40. Rossi M, Lattanzi A, Cortese L, Amodio D. An approximated computational method for fast stress reconstruction in large strain plasticity. Int J Numer Meth Eng. 2020;121(14):3048-65. https://doi.org/10.1002/nme.6346.

41. Harris CR, Jarrod Millman K, van der Walt SJ, Gommers R, Virtanen P, Cournapeau D, et al. Array programming with NumPy. Nature. 2020;585(7825):357-62. https://doi.org/10.1038/s41586-020-2649-2.

42. Nelder JA, Mead R. A simplex method for function minimization. Comput J. 1965;7(4):308-13. https://doi.org/10.1093/comjnl/7.4.308.

43. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: fundamental algorithms for scientific computing in python. Nat Methods. 2020;17:261-72. https://doi.org/10.1038/s41592-019-0686-2.

44. Takizawa H, Oide K, Suzuki K, Yamanashi T, Inoue T, Ida T, et al. Development of the user subroutine library Unified Material Model Driver for Plasticity (UMMDp) for various anisotropic yield functions. J Phys: Conf Ser. 2018;1063: 012099. https://doi.org/10.1088/1742-6596/1063/1/012099.

45. NumPy: F2PY user guide and reference manual — NumPy v1.23 Manual.

46. Martins J, Thuillier S, Andrade-Campos A. Identification of material parameters for plasticity models: a comparative study on the finite element model updating and the virtual fields method. AIP Conf Proc. 2018;1960: 110007. https://doi.org/10.1063/1.5034964.

47. Dassault Systèmes: Abaqus 2019 User's Manual.

48. Fu J, Xie W, Zhou J, Qi L. A method for the simultaneous identification of anisotropic yield and hardening constitutive parameters for sheet metal forming. Int J Mech Sci. 2020;181: 105756. https://doi.org/10.1016/j.ijmecsci.2020.105756.

49. Martins J, Thuillier S, Andrade-Campos A. Calibration of a modified Johnson-Cook model using the virtual fields method and a heterogeneous thermo-mechanical tensile test. Int J Mech Sci. 2021;202-203: 106511. https://doi.org/10.1016/j.ijmecsci.2021.106511.

50. Pierron F, Avril S, Tran VT. Extension of the virtual fields method to elasto-plastic material identification with cyclic loads and kinematic hardening. Int J Solids Struct. 2010;47(22):2993-3010. https://doi.org/10.1016/j.ijsolstr.2010.06.022.

51. Marek A, Davis FM, Pierron F. Sensitivity-based virtual fields for the non-linear virtual fields method. Comput Mech. 2017;60(3):409-31. https://doi.org/10.1007/s00466-017-1411-6.

52. Martins JMP, Thuillier S, Andrade-Campos A. Calibration of anisotropic plasticity models with an optimized heterogeneous test and the virtual fields method. In: Baldi A, Kramer SLB, Pierron F, Considine J, Bossuyt S, Hoefnagels J, editors. Residual Stress, Thermomechanics & Infrared Imaging and Inverse Problems. vol. 6 of Conference Proceedings of the Society for Experimental Mechanics Series. Springer; 2020:25-32.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.