Experimental Validation of the Sensitivity-Based Virtual Fields for Identification of Anisotropic Plasticity Models

In this work, the sensitivity-based virtual fields have been applied to identify two anisotropic plasticity models (Hill48, Yld2000-2D) using a deep-notched tensile test performed on flat samples of cold-rolled sheet of DC04 steel. The material was characterised using the standard protocol to obtain the reference sets of parameters. Deformation data was obtained during deep-notched tests using stereo digital image correlation and the virtual fields method was employed to identify material parameters. It was found that the sensitivity-based virtual fields outperform the standard user-defined virtual fields in terms of accuracy.


Introduction
To describe material behaviour accurately, models have become increasingly complex and involve more and more material parameters that need to be identified from mechanical tests. Typically, parameters are measured with a number of simple homogeneous tests, where each test provides limited information about the inferred model. As a result, many tests are generally needed to fully characterise such material models. On the other hand, developments in full-field measurements offer the ability to collect large amounts of data with the potential to improve identification of material properties. This can be used to design a new class of tests, where deformation is heterogeneous, leading to a range of multi-axial stress states within a single specimen. Probing material behaviour under such loading provides an opportunity for a reduction of the number of tests needed for characterisation, and the development of better models. The problem of increasing amount of experimental effort needed to characterise a material is an important one for the sheet metal forming community, where accurate characterisation of plastic anisotropy is essential. For instance, the simplest anisotropic model, Hill48 [1], requires three uniaxial tests performed at three distinct orientations (0 • /45 • /90 • ) to be fully characterised for plane stress applications. The model however is well known to be performing poorly, especially under biaxial loading where it cannot accurately represent the behaviour of most commonly used alloys. Many improvements have been proposed, often involving biaxial data in the formulation, significantly increasing the experimental effort involved in identifying the models. Popular models that can accurately capture the response of sheet metals, such as: Yld89 [2], Stoughton's model of 2002 [3], BBC2000 [4], BBC2005 [5] and Yld2000-2D [6] require four tests in total: three uniaxial and one equibiaxial tests. Furthermore, there are even more complex material models such as CB2001 [7,8] involving five uniaxial and one biaxial tests and Stoughton's model of 2009 [9] or Yld2004-18p [10] that require seven uniaxial (performed in increments of 15 • ) and one biaxial tests. Although these models improve the accuracy of numerical simulations, they often involve an extensive experimental effort to characterise the material.
A possible strategy to reduce the number of experiments to identify a given model is to use full-field measurements and design heterogeneous tests from which more data could be collected, compared to homogeneous counterparts.
Standard tests produce uniform/simple stress fields that can be analytically linked to the applied load, these are generally referred to as 'statically determinate'. This approach produces a single data point on the yield surface per test. Heterogeneous, statically non-determinate tests on the other hand produce a cloud of points in the stress space, each exhibiting a unique combination of stress/strain states, to which the model can be matched. In this approach, test design and material parameters extraction from the collected data represent specific challenges, as the stress field is not known a priori. To identify material parameters from full-field data, inverse techniques need to be employed.
Two of the most popular inverse techniques for extracting constitutive parameters from full-field measurements are finite element model updating (FEMU) and the virtual fields method (VFM). In the former, a model of the experiment is built up using finite element method (FEM) and the experimental data are matched to their simulated counterparts. The matching can be done based on the loading force, displacements, strains, or even the identified biaxial stress [11]. In the VFM, the stress equilibrium is enforced over the entire region of interest (ROI). It depends on the stress field reconstructed from the measured deformation through the assumed constitutive law and the material parameters are found such that they minimise the gap in the stress equilibrium. The method has successfully been applied to metal plasticity [12][13][14][15], composites [16], concrete [17], elastomers [18][19][20][21] and biological tissues [22,23] to cite but a few. One of the advantages of the VFM over FEMU is its computational efficiency; it was reported that it the VFM is 125 times faster when applied to anisotropic hyperelasticity [23] and approximately 300 times faster in the context of metal plasticity [24]. Recently a new technique emerged, called integrated digital image correlation (IDIC) which in essence combines the steps of DIC and FEMU into a single procedure to improve their metrological performance. It was applied to identification of plastic parameters by Ruybalid et al. [25], Mathieu et al. [26] and Bertin et al. [27,28].
One of challenges in performing a successful statistically indeterminate test is to ensure that it contains enough heterogeneity, i.e. a sufficiently large number of unique stress states that describe the entire constitutive model. In practice, the test is usually performed on a standard test machine and the heterogeneity is achieved by means of the geometry of the specimen. Notched samples were proven to be particularly popular for testing ductile materials. With sufficiently deep notches it is possible to activate all stress components, which is important when dealing with anisotropic materials [11,15,29]. An alternative is to machine a special specimen such as Σ-shaped sample in [14]. A methodical approach to design adequate heterogeneous tests is still an open problem. Recently, there were a few attempts at using optimisation techniques that iterate through a number of design variables to improve a measure for strain heterogeneity [30,31]. For other constitutive models, test design optimization has been studied in more depth, initially using strain heterogeneity metrics as well [32], then using balanced identification uncertainty over the whole set of parameters [33,34]. However, all failed to take into account the systematic error arising from the finite spatial resolution of the camera. The next generation of test optimization procedures relies on synthetic image deformation and minimizes the maximum identification error including the systematic error [35,36]. Extending this strategy to elastoplasticity models is the next step.
In terms of anisotropic plasticity, a number of different test configurations were used to identify popular models, with most of the effort dedicated to Hill48, due to its popularity and simplicity. The problem was tackled as early as in 1998 by the pioneering work of Meuwissen et al. [37] who used a specimen with asymmetrically placed notches. They measured displacements using a discrete number of trackers, and compared them with a numerical model to fit the parameters. Since then, many approaches have been adopted to characterise Hill48 [14,15,[27][28][29][38][39][40][41][42][43], the Ferron model [44], Yld96 [38], Bron and Besson model [42,43] and Yld2000-2D [11,15]. They included a mixture of tensile tests performed on specimens including geometric features such as holes or notches and non-standard biaxial tests leading to a heterogeneous state of stress. However, Hill48 has proved inadequate to accurately describe the behaviour of many anisotropic elastoplastic materials [9,15,42]. One of the challenges of fully characterising more advanced constitutive models such as Yld2000-2D is the activation and identification of all material parameters involved (eight or more).
In the VFM the effectiveness of parameters extraction relies on using robust virtual fields (VFs). These are spatial weight functions allowing to probe parts of the specimens for information. Traditionally, they are defined manually by the investigator using analytical functions such as polynomials, sinusoids, exponential functions etc. and are called user-defined virtual fields (UDVFs). This approach was successfully applied to the case of anisotropic plasticity by [14,15,45]. However, it was noted that the choice of virtual fields was essential for good accuracy. The choice is dependent on the expertise of the user, and might be time consuming as it involves a trial-and-error procedure. Moreover, this intuition-led choice has no reason to be optimal for the extraction of all parameters. This is particularly important for less influential parameters which may only affect the deformation fields over certain time steps and specific areas of the specimen. Recently, a new type of virtual field has been proposed to address the limitations described above, namely the sensitivity-based virtual fields (SBVFs) [46,47]. These fields are automatically generated for any constitutive model, and any test geometry based on the sensitivity of the reconstructed stress field to each material parameter. This framework provides enhanced flexibility and allows to tackle complex constitutive models more effectively [48].
In this work, we present an experimental validation of the sensitivity-based virtual fields for anisotropic plasticity. We have tested a cold-rolled sheet of DC04 steel and performed standard characterisation to obtain material parameters for Hill48 and Yld2000-2D models. Then, heterogeneous tests (deep-notched specimens) were performed along different orientations and the VFM with the SBVFs were used to characterise the models.

Brief Recall of the Finite Deformation Framework
Let us consider a body B, where the position of particles in the reference configuration is given by X and in the deformed one by x. The displacement field is defined as the difference between the current and the reference positions: The deformation gradient is defined as: where I is the second order identity tensor. Using polar decomposition, the deformation gradient can be written as the product of two second order tensors: where V is the left stretch tensor and R is the rotation tensor. The left stretch tensor can be conveniently calculated as: where the root operator refers to the root of a matrix.
A consequence of such mathematical description is that for every point, a local coordinate system rotates during deformation, as outlined in Fig. 1. This is an important feature to consider when the body includes a texture, as its orientation will follow any local rotations. A convenient measure of strain, called Hencky strain, can be constructed from the left stretch tensor: This strain measure can be used to formulate constitutive laws within the finite deformation framework. For further details on continuum mechanics the reader is referred to [49]. Definition of coordinate systems, X is the initial position of a material point and x, its current position, (i, j ) is the initial orientation of the local coordinate system, (1, 2) is the corotational system, (Ξ, H ) is the material coordinate system in the reference configuration and (ξ, η) is the material coordinate system in the current configuration

Constitutive Models
In this study we considered two different yield models suitable for cold-rolled sheets: Hill48 and Yld2000-2D [1,6]. The former is relatively simple extension of the von Mises criterion to account for anisotropy and the equivalent stress can be expressed as: The coefficients in this criterion can be defined in multiple different ways. Here, we follow the convention used in finite element package Abaqus. For alternative ways of defining the coefficients the reader can refer to the appendix. The plastic potentials (R ij ) are generally used to obtain the governing parameters from an experiment [50]: where R ij = σ y ij σ 0 , with σ y 11 and σ y 22 the yield stresses identified in planar uniaxial tests conducted at 0 • and 90 • respectively and σ y 33 the through-thickness yield stress. Finally, σ y 12 is the yield stress identified under pure shear. Although the model is defined for plane stress and σ 33 = 0, the information about σ y 33 can be obtained from the combination of associated flow rule assumption and Lankford coefficients. The reference yield stress was assumed to coincide with σ y 11 , i.e. σ 0 = σ y 11 , as this reduces the number of variables to be identified by one, and does not affect the formulation of the model. Additionally, plane stress was assumed, as the tested samples were thin relative to their in-plane dimensions, and associated flow rule was used.
Yld2000-2D was developed strictly for plane stress condition for which the equivalent stress can be calculated as: where a is an exponent based on the metal micro-structure (a = 8 for FCC and a = 6 for BCC) and X 1 , X 2 and X 1 , X 2 are the principal values of two stress tensors X , X which are defined as linear combinations of the Cauchy stress: Matrices L and L are given by: The model involves 8 independent parameters, α 1 -α 8 , which are generally obtained using 3 tensile tests performed at 0 • , 45 • , 90 • , a biaxial test (bulge test), as well as a test for measuring r b =ε yẏ ε xx at balanced equibiaxial loading. In total, eight parameters have to be measured to calibrate α parameters; four yield stresses: σ 0 , σ 45 , σ 90 , σ b , and four Lankford coefficients (r) characterising anisotropy in plastic deformation: r 0 , r 45 , r 90 , r b . Finally, the associated flow rule was assumed here as well.
A non-linear isotropic hardening power law (Ludwik) was chosen with the following form: where σ 0 , K, n are material parameters andε p is the equivalent plastic strain integrated throughout the history of deformation.
Reconstruction of the stress tensor from experimentally measured strain data is performed using a numerical implementation of the constitutive laws. These were based on the radial-return algorithms developed by Koh et al. for plane stress Hoffman model (generalisation of Hill48) and Yoon et al. for Yld2000-2D [51,52]. Note that the Hoffman model simplifies to Hill48 when tensioncompression symmetry is assumed.
The constitutive computations are performed in a material coordinate system (ξ, η) which is initially aligned with manufacturing rolling (RD) and transverse (TD) directions. Since the kinematic fields are computed in the global frame (i, j ), for each data point strain and stress tensors need to be rotated to the material frame: where R mat is a rotation tensor projecting the global frame onto the material frame in the unloaded configuration ( Fig. 1). Once the stress tensor is reconstructed it is rotated back to the global frame in which the VFM equations are formulated.

Virtual Fields Method
The Virtual Fields Method is an inverse technique to identify material parameters from full-field measurements. It relies on the force equilibrium through the principle of virtual work (PVW). In the case of static loading and in absence of body forces, it can be expressed in the reference body configuration as [53]: where B 0 is the considered body in the reference configuration, ∂B 0 its boundary, N is the outward vector of ∂B 0 , P is the 1 st Piola-Kirchhoff stress tensor and U * is a vectorial test function called virtual displacement. Virtual displacement fields need to be continuous and piecewise differentiable. The stress tensor (P) is directly reconstructed from measured kinematic fields by means of the assumed constitutive law and a guessed set of material parameters (χ). The validity of the guess is assessed by the residual value evaluated with equation (14).
As the full-field measurements provide spatially rich data, the first integral in equation (14) can be replaced by a discrete sum of all points in the region of interest (ROI). By selecting virtual displacements as constant across ∂B 0 , the second integral in equation (14) can be replaced by the product U * · F load , where F load is the total load measured with the test machine load cell. This procedure filters out generally unknown distribution of tractions over the loading edges with a quantity easy to measure.
Since the PVW has to be evaluated over the entire volume, but the measurements are taken only at the surface of the specimen, it is necessary to make some assumption regarding the variation of kinematic fields through the thickness. If the specimen is thin, plane stress assumption is reasonable and the PVW can be expressed in terms of a cost function as: where α i is a scaling parameter to normalise the contribution of each virtual field and their respective units, S j is the surface area of each measurement point and h is the thickness of specimen. This cost function can include a number of independent virtual fields and load levels (time steps).
The identification is carried out by minimising (equation (15)) with respect to the sought material parameters. It is worth noting that equation (14) is formulated in terms of the 1 st Piola-Kirchhoff stress tensor, while most constitutive laws relate kinematic fields to the Cauchy stress tensor. The former can be obtained from the latter with: where F is the deformation gradient tensor. Since in reality the deformation is fully 3D, so is the deformation gradient and this has to be reflected in the computation of its determinant. By assuming negligible out-of-plane shearing, the determinant can be computed as: Since the in-plane values are directly measured the only unknown is the out-of-plane component. It can be estimated from the constitutive law (e.g. assuming plane stress elasticity and isochoric plastic flow) [47], or can be directly measured by means of back-to-back camera systems as shown in [54]. The virtual displacements (and their spatial derivatives, later simply referred to as virtual fields) spatially probe the reconstructed stress field for information. As a result, the choice of VFs is crucial and has strong influence on identification quality. Ideally, VFs should be focused on areas rich in information and minimize the influence of the measurement noise. In the case of non-linear material models, selecting VFs usually relies on a manual definition by means of simple mathematical functions, such as polynomials or sinusoids [14,15]. The effectiveness of these VFs called user-defined VFs (UDVFs) heavily depends on the expertise of the investigator. It is worth noting that usually UDVFs are kept constant across the history of deformation, whereas the information evolves as the loading changes. As a result, a VF that is relevant e.g. for identification of yield-related parameters might not be as effective for identification of the hardening law. Recently, a new automated method for generating highquality virtual fields has been proposed [46,47]. These fields called sensitivity-based virtual fields (SBVFs), are designed to highlight areas rich in information for each parameter separately and for each time step, resulting in better identification, without significant input from the investigator.

Sensitivity-Based Virtual Fields
The sensitivity-based virtual fields [46,47] are automatically generated for every material parameter. These fields are good at finding information about each parameter separately and coupling them together in one cost function to identify all material parameters. Each virtual field is constructed based on the sensitivity of the reconstructed stress field to a given material parameter. The procedure for generating virtual fields is discussed in details in [46] with the extension to large deformation framework in [47]. Here, we shortly summarize the necessary steps to generate SBVFs.
For each material parameter a map of stress sensitivity is calculated through: where δχ i is a small variation of parameter χ i , typically δχ i is between -20% and -10% of χ i . The fields are calculated at all considered load levels, generating temporal maps of stress sensitivities. In plasticity, the response is history dependent with yieldrelated parameters being active from the onset of plasticity throughout the whole history, while hardening parameters are active only during accumulation of plastic deformation. In order to decouple the influence of yield stress from that of hardening, an alternative stress sensitivity field is derived, called the incremental stress sensitivity field: These fields highlight the information about a given parameter in the test and are excellent candidates for virtual fields. In order to apply them in the VFM, virtual fields are generated such that the spatial derivatives of virtual displacement fields match the incremental stress sensitivity fields in a least-square sense. Additionally, the VFs need to be constructed in a way that the corresponding virtual displacement field is known. To achieve that, a virtual mesh is employed to express virtual fields using piecewise linear functions. It consists of linear quadrilateral elements enclosing the ROI, which are used to express both virtual displacement and virtual strains fields based on the nodal values through interpolation functions.
To construct SBVFs, a virtual global strain-displacement map is constructed: for each data point, the spatial derivatives of the virtual displacement at that point are written as functions of the nodal virtual displacements through the element shape functions. All such equations are concatenated in one matrix (B glob ). The matrix is then modified to account for the wanted virtual boundary conditions (e.g. U * = const over ∂B 0 ), yielding the modified virtual strain-displacement matrix:B glob . Then, the incremental stress sensitivity map is projected in a leastsquare sense onto the virtual mesh, using a pseudo-inverse ofB glob : Finally, the identified virtual displacements are used to calculate virtual fields at all data points: The procedure yields a separate virtual field for every material parameter, at every data point and every time step which is then used to evaluate (equation (15)). Equation (20) leads to virtual displacements that have a unit of stress sensitivities multiplied by length. However, it needs to be reminded that the virtual displacements have no physical significance and the unit is unimportant. For a better illustration, one can imagine a separate field is constructed, which is identical to the stress sensitivity field in terms of numerical values but unit-less. This field is then projected onto the virtual mesh leading to virtual strains with the unit of 1/length, which is consistent with the principle of virtual work. When multiple virtual fields are combined in one cost function, virtual displacements need to be normalised so that they have comparable weights. This can be done by tuning the α i parameter in equation (15) in a way to normalise the magnitude of e.g. the internal virtual work. This was done previously based on the mean value of x% mostly contributing time steps as described in [46,47].

Specimen Preparation
Specimens were water-jet cut out of a cold-rolled sheet of DC04 low-carbon steel alloy with a nominal thickness of 1.5 mm. Three geometries were tested: a standard dogbone (DB) specimen, a rectangular specimen for bulge test and a deep-notched (DN) specimen for heterogeneous test (Fig. 2). The specimens surfaces were first cleaned with sandpaper to remove oxides and then coated with a rubberbased white paint (Rust-oleum Peel coat, white matt finish) to provide good contrast for black speckles. An optimised speckle pattern [55] was printed on the specimen using a flat-bed printer (Canon Océ Arizona 1260 XT) which provided good consistency and reproducibility. An average speckle size of approximately 65 μm was achieved (Fig. 3).

Tensile testing
Both DB and DN samples were tested using a servohydraulic test machine with a 100 kN load cell and hydraulic grips. Deformation was measured using a stereo-digital image correlation (DIC) set-up, with two digital Manta G-504b cameras (5 Mpx), equipped with 105 mm Sigma DG Macro lenses and polarisers (Edmund optics). A LED light panel (Zaila, Nila) equipped with a polariser was used to illuminate the samples. By setting cross-polarization specular reflection was minimised which resulted in a greylevel histogram spread across most of the dynamic range of the cameras [56]. The stereo-DIC setup was used following the recommendation of the DIC Good Practices Guide, to account for test piece thinning, rotation or translations induced by misaligned grips and other effects that may produce errors in 2D-DIC [57]. The reference images for correlation were taken while maintaining zero load; the specimen was loaded in displacement control divided into three phases: slow rate (elastic range), medium rate (transition to high rate) and high rate where most of the plastic deformation took place, as indicated in Fig. 4. The length and rate of displacement of each phase was tuned in the preliminary phase of the testing campaign and were proven to give a good compromise between the extent of collected data and the duration of the test. The images were taken every 1 s and synchronised with the force measured   Table 1.

Bulge tests
Since the yield stress and r-value information in balanced biaxial state is required for determination of the anisotropic constitutive parameters, a hydraulic bulge test was carried out using an Erichsen bulge/FLC tester model 161.
A steel sheet specimen was mounted on the bulge test apparatus as in Fig. 6. Then the specimen was clamped between the lower blank holder and the upper die. To prevent slipping of the specimen during the test, plastic flow was restricted with a drawbead and high blank holding force. The diameter of the area of interest in the test device was 200 mm. The hydraulic pressure was applied on the bottom side of the specimen to produce bulging and plastic deformation.
A stereo digital image correlation technique was used for measurement of curvature and strain fields as shown in Fig. 7. Two 2448 × 2048 pixels 14 bit CCD cameras were used for the measurement. The important correlation variables chosen in the DIC analysis were: subset: 41, step: 7, object pixel size: 0.117 mm; for the details the reader is referred to the DIC table in Appendix C.  In the bulge test, the biaxial stress-strain curve is derived from the membrane stress and the through-thickness strain near the pole of the bulged specimen according to the procedure outlined in [58].
In membrane theory of a thin-walled pressure vessel, the biaxial stress is calculated as: where p is the pressure obtained from a pressure sensor, R the current radius of curvature and t the current thickness. The current thickness is calculated from: where t 0 is the initial thickness and ε t the through-thickness strain, obtained from the in-plane strains through: ε t = −ε 1 − ε 2 , where ε 1 and ε 2 are major and minor strains. The current radius of curvature R is obtained from: where κ is the curvature. Both the in-plane strains and curvature were measured with Vic-3D (Correlated Solutions). The in-plane strains were calculated by averaging major and  minor strains within the circle centred at the apex and with the radius r 2 = 10 mm. The curvature κ was calculated by averaging curvatures κ x and κ y along the horizontal and vertical diameters of a circle centred at the apex and with the radius r 1 = 20 mm respectively. Two specimens were tested and it was found that the deviation was very small between the two membrane stressthickness strain curves.

Dogbone specimens
Raw grey-level images were exported to a DIC package (MatchID 2018.2.2) and processed using stereo-DIC. Due to significant plastic deformation of the specimen upstream of the gauge section the small ROI that remained in the camera field of view for the whole test was selected. The camera parameters and DIC settings are summarized in Table 8 available in Appendix, as per recommendation from the International DIC Society Good practice guide [57]. The measured fields were used to reconstruct the stress-strain curve and identify yield stress and hardening law for each of orientation (0 • /45 • /90 • ), using the uniform and uniaxial stress assumptions. The average longitudinal plastic strain was plotted against the average transverse plastic strain, and a straight line was fitted to the data in order to obtain the Lankford coefficient [59]. The line passed through the origin and was fitted to the data corresponding to 8-12% range of plastic deformation.

Deep-notched specimen
Raw grey-level images were exported to MatchID and displacements in the gauge section were obtained using stereo-DIC, with the parameters presented in Table 9 included in Appendix C. The data was exported to Matlab, where displacements were temporarily and spatially smoothed and then down-sampled to a number of load levels (time steps). Gaussian filter (defined with the standard deviation σ smooth and the kernel size is chosen as the lowest odd number that is larger than 6 × σ smooth ) with edge correction was used for spatial smoothing and Savitzky-Golay filter was used for temporal smoothing, characterised with two parameters: polynomial order m T S and the window size w T S . The down-sampling purpose was to improve signal-to-noise ratio of strain increments and reduce the computational Fig. 9 Longitudinal versus transverse plastic strains measured using dogbone samples. The data was used to calculate Lankford coefficients in the range of 8%-12% effort. Central finite difference was used to compute kinematic data (deformation gradient, rotation tensor, Hencky logarithmic strains), which were then passed to an in-house VFM code to identify the material parameters.

Dogbone Testing
For each of the orientations (0 • /45 • /90 • ), three specimens were tested. The measured true stress-strain curves are presented in Fig. 8. A visible bump around the strain value of 3% is believed to be caused by the step change in the cross-head velocity as indicated in Fig. 5 and is consistently seen in all collected data (including DN samples). The curves were used to identify the hardening parameters (equation (12)) and coefficients for Hill48 (equation (6)) and Yld2000-2D (equation (8)). The hardening law was identified using strain of up to 10%, as at the higher deformation the assumed hardening model does not capture the material behaviour accurately. This could be improved upon in the future, however the main objective of this contribution is to demonstrate the effectiveness of the SBVFs, as opposed to improving the constitutive description of the material. Apart from the yield stress defining the hardening  Longitudinal plastic strain was plotted against transverse plastic strain as shown in Fig. 9. The Lankford coefficients were identified in the range of 8%-12%. Across all samples the trend has shown good linear relationship with little variation between different samples. All of the identified parameters are presented in Table 2.
The reference hardening curve was constructed based on the average behaviour along the RD: The ratios of flow stress in each direction to the flow stress in RD have been investigated under different values of plastic work. This is now a standard practice in evaluating anisotropy of metal sheets, as the ratios are difficult to identify reliably at low levels of plastic deformation, especially when the bulge test is employed [10]. The obtained curves are presented in Fig. 10. 1 Clearly, the anisotropy in yield stress changes rapidly at low deformation levels. This has to be taken into consideration when selecting the reference parameters for comparison with the DN results. It is reasonable to consider the average ratios based on plastic work between 10 and The collected data was used to obtain Hill48 parameters in two ways: using all Lankford coefficients, or using σ 0 from all tests and the Lankford coefficient from DB-90, for the details the reader is referred to the appendix. The variation of the initial yield stress and Lankford coefficient with material orientation is presented in Fig. 11. The corresponding parameters are presented in Table 3. The two sets are different due to the limitations of the model and its inability to describe different degrees of anisotropy of the yield stress and deformation at the same time [60]. Importantly, the two sets differ only by the magnitude of σ y 12 which drives the values at 45 • .
The biaxial flow stress (σ b ) and the strain ratio at balanced biaxial loading (r b = 0.77) were measured in the bulge test. The ratio of σ b to the RD flow stress as a function of work hardening is presented in Fig. 10. The equibiaxial yield stress varies significantly at low level of deformation due to differential hardening [61]. This ultimately leads to a set of different Yld2000-2D surfaces identified from the standard tests, two of the surfaces identified (with ratios taken at 10 and 25 MPa) are presented in Fig. 12. The method used for fitting the parameters is detailed in the appendix. As seen from the plot, the biaxial yield stress is over-predicted when Hill48 is used, due to lack of flexibility of the model. Naturally, Yld2000-2D offers a much better fit to the experimental data as demonstrated in Fig. 11, where the variation in both yield stress and Lankford parameters is captured correctly with a single set of parameters.
The identified parameters for the Yld2000-2D models are presented in Table 4 for two different levels of work hardening (10 MPa and 25 MPa). Exp Mech (2020) : -

Fig. 12
Identified yield surfaces using the standard testing protocol using two models: Hill48 and Yld2000-2D

Deep-Notched Test
Three directions were tested: (30 • /45 • /60 • ), so that different test configurations could be combined together in a single cost function. For instance, by combining two tests performed at 30 • and 60 • , richer data is available to identify the models, which should improve accuracy of identification. Additionally, these other orientations may be used as a validation tool for testing the performance of the model and the quality of identified parameters. As indicated in Tables 8-9, DN tests exhibited a smaller strain uncertainty. This is believed to be caused by a better lighting set-up during DN experiments. Figure 13a 2 shows an image of a typical specimen in the unloaded configuration. Only the region bounded by the solid box was correlated and the data corresponding to the dashed box was then passed on to the identification routine. The magnification in the experiment was chosen such that the field of view used in identification would be captured by both of the cameras, regardless of the deformation of the ROI or the rigid body motion due to the deformation upstream the test piece. Additionally, if the cameras were positioned closer to the test machine, obtaining sufficient depth of field would be challenging. The strain fields obtained at 10.5 kN (45 • specimen) are presented in Figs. 13b-d and the force-displacement curves 2 The spikes in the histogram were identified after the experiments were conducted and come from a (now resolved) software issue when saving 8-bit images. Since the quality indicators (grey level noise, displacement and strain noise-floors) seem to be reasonable it is suspected that those artefacts do not significantly affect the measurements. are shown in Fig. 14. The displacements in that plot were first corrected for rigid body motion and then measured at the point indicated in Fig. 13(a) with a red dot. The figure demonstrates that consistent measurements were obtained with multiple test specimens and that the variation in force measurements between the three orientation is minor. Note, that although force-displacement curves for the three orientations are similar, they correspond to different deformation fields.
Samples were strained until failure and the base of white paint did not fail during the experiment, however it debonded from the specimen when the neck started to develop, see Fig. 15. If observation of the strain localization was of interest, the base layer could be removed, and white speckles could be used instead of black ones. In that case, the contrast could be achieved by means of the material surface, combined with cross-polarization of light to remove the specular reflection.
Displacements obtained in MatchID were exported to Matlab, temporarily smoothed with m T S = 3 and window w T S = 11, then the data was down-sampled to the desired number of load levels and spatially smoothed with a Gaussian of square kernel of 13 × 13 pixels (corresponding to σ smooth = 2).
A representation of the data in minor-major strain space is shown in Fig. 16. DN tests occupy the space around uniaxial tension, with a large number of points lying in between uniaxial tension and pure shear/plane strain. As expected, the data is contained in the second quadrant and thus not representative of biaxial tension which is relevant to the Yld2000-2D model.
The kinematic data was fed to the in-house VFM program to identify the parameters of Hill48 and Yld2000-2D. For the identification, a virtual mesh of 10 × 10 elements, a material variation of 15% (δχ i = −0.15χ i ) and a scaling parameter (see [47] for more details) of 0.3 were used. As the baseline total stress sensitivities were used to construct SBVFs as they yielded better results, as shown later. Additionally, the elastic properties were set a priori to the typical values for steel, i.e. E =200 GPa and ν = 0.3.
Minimisation of the cost function (equation (15)) was carried out in Matlab, starting from four points selected with a random number generator, from the region contained between 50% and 200% of the reference values. Initially, the Levenberg-Marquardt algorithm was used (lsqnonlin) but it was found that the algorithm could not converge to a unique global minimum. The method was then switched to fmincon with the sequential quadratic programming algorithm (SQP) which was capable of converging to the same solution regardless of the starting point. The gradient of the cost function was calculated internally by Matlab using the central finite difference scheme. To quantify the identification accuracy a metric based on the reconstructed apparent yield stress was used [47]. In this approach, instead of comparing parameters on a one-to-one basis, the models are used to reconstruct the apparent yield stress that would cause yielding at a given orientation in a 1D loading scenario (see Fig. 11a). Additionally, the effect of hardening can be accounted for which results in a map of the apparent yield stress as a function of orientation and level of plastic strain. The global mean error is constructed as a mean difference between the map corresponding to parameters derived from the dogbone tests and the map corresponding to the identified parameters. This procedure helps quantify the identification error in a more meaningful way than the standard one-to-one parameter comparison.

Identification of Hill48 with SBVFs
The first set of parameters was identified using specimens cut at 45 • to RD. In total, four different samples were tested and processed, but the data from DN-45-2 was discarded as it was found that there was significant out-ofplane bending during the test which violates the throughthickness homogeneity assumption. The data from the three successful tests were passed to the in-house VFM code and the details regarding time steps included in the identifications are presented in Table 7. Note that frames were selected such that the majority of points were within 10% of strain in the load direction so that the DN data were consistent with the hardening law definition assumed for the  dogbone tests. The criterion for the final time step was that the strain corresponding to μ + 2σ was about 10%, with μ being the mean and σ the standard deviation of the strain in the loading direction. This ensures that about 98% of data is below 10% strain. On average each identification run took 3 hours to complete. The identified parameters were used to compute the mean error based on the apparent yield stress metric. A typical map corresponding to DN-45-3 is shown in Fig. 17. The dashed cut (Fig. 17b) shows how well the yield stress was reconstructed at equivalent plastic strain of 0.002, and the solid cut shows the hardening behaviour along the orientation of 45 • , indicating good identification of both anisotropy and hardening. The identified parameters and the corresponding mean error values for all samples are shown in Table 5, and graphically in terms of the apparent yield stress and predicted Lankford coefficients in Fig. 18.
The results show that the parameters obtained with the DN tests are repeatable and in good agreement with the dogbone data in terms of yield stress and hardening, however the Lankford coefficients are significantly lower than for the dogbone.The vertical position of the reconstructed  Lankford curves is driven by σ y 33 which also controls how far the biaxial part of the yield surface extends. As in the DN tests much richer information is available in comparison the homogeneous counterparts, the limitations of the model become much clearer. Because of the interactions between the yield surface, the flow potential and the biaxial yield stress, it is impossible to match all the data at the same time with Hill48. This signifies that more advanced constitutive models are required to accurately describe the material under investigation.

Identification of Hill48 with UDVFs
An alternative to the SBVFs are the user-defined VFs. A set of viable virtual fields for Hill48 and the geometry used in this test has been presented in [15] based on a trial-anderror method and expertise of the lead author. The virtual displacements were defined in the coordinate system presented in Fig. 19 and were constructed in a way to include all stress components in the cost function. These fields are: The same data as used in Section "Identification of Hill48 with SBVFs" were fed to the VFM algorithm, however now the UDVFs were used instead. The identified parameters were used to reconstruct the variations of the yield stress and Lankford coefficient with orientation as shown in Fig. 20.    It is apparent from the figure that the UDVFs were not as successful at characterising Hill48 model as SBVFs were. The yield stress was accurately predicted at about 45 • which corresponds to the orientation of the test specimen, however large differences were noted elsewhere. By using the same metric as before, the global difference was calculated and was found to be in the range 6%-19%. Judging on the Lankford coefficients it is clear that the reconstructed parameters are not physical as they lead to a non-continuous distribution.
The first virtual field (equation (26)) represents a uniform extension, which leads to a direct comparison between the measured force and the force reconstructed from the stress field. In the case of anisotropic properties this integral quantity is not sufficient to describe the model. Although other virtual fields have been included in the cost function, they were formulated in a way not to include the work of external forces. As a result, the residuals corresponding to those fields are much smaller than the one of the first field and they are not very effective at including the other two stress components in the cost function. These results highlight the problem of manually defining virtual fields, as they need to be hand tailored to every application with great care and expertise in order to be functional.
The parameters identified in DN-45 tests with both UDVFs and SBVFs were used to predict the internal reaction forces based on the stress field reconstructed from the measured strains. These forces were compared to the experimentally measured ones to give an independent validation in terms of quality of the material model and the accuracy of identification with the DN-30 and DN-60 tests. Fig. 21 shows the reconstructed forces corresponding to the properties identified with DN-45-4 and strain measurements from DN-45-4, DN-30-1 and DN-60-1. In the case of DN-45-4 both SBVFs and UDVFs accurately reconstructed the force measurement. When these parameters were used to predict the responses of DN-30-1 and DN-60-1 a large discrepancy was found with UDVFs, indicating poor identification of anisotropy by UDVFs. On the contrary, the parameters derived from the same test, but using SBVFs, lead to a much better prediction of the experimentally measured force-displacement curves indicating that a single test performed along 45 • is sufficient to identify the material plastic anisotropy with the Hill48 model, which confirms the results obtained in [27].

Influence of data range on Hill48 parameters identification
In order to investigate whether the results depend on how much data is used in the identification, a study was performed where increasingly more time steps were fed to the cost function. Specimen DN-45-4 was chosen as the baseline, which in the original study contained the maximal strain of 14.0%.
Here, we investigated how the identified parameters depend on the number of frames used in the identification. At the low end, the maximum strain in the load direction was 9.6%, and at the high end 20.3%, with the intermediate steps approximately every 1.4% of additional strain, as indicated in Fig. 22. In principle, as larger deformations are supplied to the identification, more emphasis should be put on the model to match the Lankford coefficients more accurately. However, the incremental stress sensitivities related to yielding tend to filter out time steps at which most of the deformation takes place. To address that, two studies were run: one in which incremental stress sensitivities were used as a base for generating SBVFs, and another one where total stress sensitivities (equation (18)) were used instead. The findings of this study are presented in Figs. 23 and 24 for incremental and total stress sensitivities respectively and the mean error as a function of maximal strain supplied is presented in Fig. 25. Surprisingly, the results obtained using the total (as opposed to incremental) stress sensitivities are much more consistent and stable with respect to the total number of time steps used, as indicated by the lower mean error over the range of maximal plastic strains. This justifies using total stress sensitivity as the baseline in this contribution and comes as surprise as previous studies showed that the incremental sensitivities generally yielded more accurate parameters. These were however performed on synthetic data so did not include any modelling errors [46,47]. Although the identification was improved when the SBVFs were calculated with total stress sensitivities, no improvement in reconstruction of Lankford coefficients was found, despite supplying total strain as large as 20%.
A possible explanation for the under-performance of the SBVFs based on incremental stress sensitivities is due to filtering out large amounts of data from the cost function by the fields. As mentioned in Section "Sensitivity-Based Virtual Fields", the incremental form of stress sensitivities was introduced to filter out history dependence of plasticity models and highlight transition zones between elasticity and plasticity. As a result, when most of the field of view (FOV) is plastically deformed, the incremental stress sensitivity fields only cover a very small part of the FOV and provide little information to the cost function. This might have a detrimental effect in anisotropic plasticity, particularly when associated-flow models are used, as the information about the anisotropy can be queried at different levels of plastic work, not only during initial yielding, but also during accumulation of deformation. This shows that more research is required into understanding of the benefits and limitations of incremental versus total stress sensitivity for the SBVFs, so that optimal choices can be made.

Identification of Yld2000-2D with SBVFs
Yld2000-2D model was identified using data from 3 experiments (DN-30, DN-45, DN-60 (25)) was used, so that only anisotropic coefficients were identified. The total stress sensitivity fields were used to generate SBVFs as there was no need to decouple yield and hardening related parameters. The cost function was minimised using the fmincon function in Matlab with the SQP algorithm. An additional constraint was added to the minimisation problem to ensure that the flow curve represented the behaviour of the material along RD [11]: Four starting points were selected using a random number generator with a lower bound of 0.5 and an upper bound of 2.0. The identification procedure was restarted after it converged for the first time to ensure that the global minimum was found. Two of the starting points converged to the same solution, shown in Table 6, and in terms of the yield surface in Fig. 26. Here, we present only results obtained with data from Set 1. The identified yield surface shows a good agreement with the reference one in the regions where data were present, however deviates significantly near the biaxial stress state. In the region covered   with data points, the difference between the two surfaces is less than 5% in terms of the shear stress. This is further confirmed by looking at the comparison between the experimental and reconstructed forces which match closely ( Fig. 26(b)). The poor description around the biaxial stress state indicates that the tests are not sufficiently heterogeneous for how flexible the Yld2000-2D function is, i.e. the data used in the identification is not diverse enough to constrain the shape of the yield surface over all possible loading states. A similar observation was recently reported in [62], where the authors demonstrated that their model was matched very well in the domain represented in the tests, however it did not provide good predictions outside of it. They suggested that adding an additional information to the cost function (in their case a test at different load rate) could significantly improve predictions over a wider domain and relieve the issue of non-uniqueness of the material parameters.
To confirm this hypothesis another identification has been carried out, using exactly the same test configurations as above, however adding an additional constraint to the optimisation problem. The optimisation function was constrained such that the identified yield surface passed through σ b and r b measured in the bulge test. Four randomly selected starting points were used, and at least two converged to a similar solution, with similar final values of the cost function ( Table 6). The minimisation procedure took between 2-3 days depending on the starting point, irrespective of whether UDVFs or SBVFs were used. The majority of the computation time was taken by the stress reconstruction algorithm, due to the large number of data points and combination of three tests. The identified yield surface, apparent yield stress and Lankford coefficients are presented in Fig. 27. The identified model matches the reference well in terms of the yield surface and the Lankford coefficients, confirming that the three deep-notched tests did not represent enough data to identify Yld2000-2D, however once the additional constraint was added to the optimisation the SBVFs did a good job at identifying the correct parameters. Importantly, both data sets were reliable as they provide consistent material parameters.
To get a further insight into the completeness of the test, the number of DN tests has been gradually reduced, first to two tests (30 • + 60 • ) and then to a single test (45 • ), both of which included information from the bulge test. The parameters were identified with these reduced tests using the same methodology as described before. To quantify the accuracy of identification, 3D yield surfaces were compared; each surface was reconstructed in the spherical coordinate system consisting of three variables: the azimuth angle, β, in σ 11 -σ 22 plane, the elevation angle, ψ, above the plane, and Table 6 Anisotropic parameters and the final value of cost function for Yld2000-2D model identified with DN tests Test configuration  26 Identification of Yld2000-2D using 3 DN tests. a Comparison of the identified Yld2000-2D yield surface and the references one (corresponding to 10 MPa of plastic work). The yield surface and the experimental data points were plotted for equivalent plastic strain of 3% to represent only the load-paths achieving significant plastic strains. The outlines are drawn for change in the shear stress corresponding to 5% of the yield stress under pure shear. b Comparison between the experimental and the reconstructed forces using identified material parameters and measured strains the radius d. For every combination of the two angles, the distance from the origin to the yield surface was calculated and compared between the reference and identified sets of parameters. We then used the relative difference between the two distances to quantify the matching: The difference metric is convenient for comparing 3D yield surfaces as it offers an one-to-one correspondence between the models. The difference maps associated with the three test combinations are presented in Fig. 28. From the figure, it is clear that the larger the number of tests included in the identification, the smaller the difference between the reference and identified models. In general, the differences were small for three and two tests set-ups (mean difference < 1%). Surprisingly, the discrepancies were most pronounced for the stress states represented in the tests. Potentially, this might be due to better fitting of the yield surface to those represented stress states, in comparison to the parameters derived from the homogeneous tests. We hypothesise that the reference set of parameters is not optimal for all stress states, since Yld2000-2D is fitted to experimental data in a least-square sense. Thus, it is possible that the deep-notch tests actually do a better job at capturing the material response under these multi-axial loading conditions. Significant differences were also found close to the uniaxial tension at 90 • . This however could be explained by the fact that the reference set of parameters was generated with the average flow stress ratios; in reality, the reference yield surface represents material behaviour approximately over a large range of plastic work, whereas DN tests contain information corresponding mostly to the low plastic work at which the ratios vary considerably (Fig. 10). In the case of the single test, the error in this region is much larger, suggesting that the test is not sufficient to characterise all possible orientation angles of the material. This was confirmed by reconstructing force-displacement curves, similarly to the procedure employed in Fig. 21. The forces corresponding to 30 • and 45 • were reconstructed with high accuracy, as they were constrained with experimental data and the optimisation constraint of equation (29). However, this was not the case for 60 • test, in which the force was predicted with an error of about 5% (data not shown).
To further elaborate on the viability of the three different test combinations, the corresponding Lankford coefficients have been calculated and are presented in Fig. 28d. Interestingly, in spite of similar mean difference levels in terms of stress envelope (Fig. 28a-b), the two tests set-up experienced larger errors in terms of Lankford coefficients close to 90 • compared to the three tests combination. This is most likely due to the fact that plastic flow is related to gradients of stress field, making it much more sensitive to the identified parameters values and experimental noise as well. It must be noted however that the information about 3D plastic deformation is present during the identification procedure as it is encoded in the strain measurements when isochoric plastic flow is assumed.
The results presented in this contribution might seem inferior to those reported by Rossi et al. [15] who used the same test to identify Hill48 and Yld2000-2D models for a stainless steel alloy. However, under closer inspection this is not the case either for Hill48 or Yld2000-2D. In the former case, they managed to correctly predict experimental forces for all three orientations only when all three tests were used to identify Hill48 parameters. As shown in Fig. 21 SBVFs achieved that successfully with only a single test. In the case of Yld2000-2D, they demonstrated that the model significantly decreases the error on Lankford coefficients reconstructions, however the error on the values they reported is of similar magnitude (Δr ≈ 0.10) to that reported in this contribution. Finally, they did not measure the reference behaviour under biaxial loading thus the quality of Yld2000-2D identification could not be fully judged. In terms of experimental forces they observed a good agreement between the predicted and measured ones. In this report we were able to reproduce this observation even with the data set not including the bulge test (see Fig. 26b), indicating that those quantities can be matched correctly while the biaxial response of the material is not matched closely. Interestingly in that paper the authors used a direct method for integrating plasticity equations which rely on the measured strains to estimate the flow direction. This procedure might process the experimental noise, modelling errors and interact with virtual fields differently leading to different results compared to the elastic predictor-plastic corrector schemes. In the future it would be interesting to investigate which of the two schemes is more suitable for identifying plastic laws.
The results suggest that the DN tests alone are not rich enough to fully identify Yld2000-2D model. It may be possible to work on the specimen shape to improve the situation but it is unlikely that this will ever be enough to place points close enough to the biaxial stress state. However, using two axes for load introduction may lead to enough degrees of freedom to obtain complete model. A test like that in [63], combined with non proportional loading may populate the stress space widely enough to obtain all parameters in one test.
Full-field measurements give an insight into a large range of load paths which the model is calibrated to. These identified parameters represent a more complete response of the material, compared to the parameters derived from three uniaxial tests and a biaxial one. Further work needs to be done to establish which of the two sets of parameters reflects the behaviour of materials more comprehensively. The validation procedure could be performed by modelling an independent test using the two models, and comparing the predictions with measurements. This is an exciting opportunity for future studies.

Identification of Yld2000-2D with UDVFs
To complete the evaluation of SBVFs performance an identification of Yld2000-2D was carried out using the UDVFs defined in Section "Identification of Hill48 with UDVFs". In this experiment the data from Set 1 (including bulge test) were used to obtain Yld2000-2D parameters. The identified parameters are shown in Table 6 and graphically in Fig. 27.
Clearly, the obtained parameters with UDVFs are not as good as the ones obtained with SBVFs, confirming that the proposed virtual fields significantly improve the identification process. The final value of the cost function is much larger compared to the SBVFs, however this was expected as the virtual fields have completely different values and scale the stress fields differently.

Conclusions
In this work we have tested a sheet of DC04 steel alloy using the standard testing protocol and experiments on deep-notched specimens to identify two constitutive models: Hill48 and Yld2000-2D. The virtual fields method, combined with the sensitivity-based virtual fields, has been employed to extract constitutive parameters from full-field measurements. The main outcomes of this study can be summarized as follows.
-The results presented in this contribution suggest that the sensitivity-based virtual fields are effective at extracting information about material parameters, particularly for advanced models such as Yld2000-2D. -When the heterogeneous tests were used to identify Hill48, the material parameters matched the yield stress variation and the loading force well, however underestimated the Lankford coefficients. We suggest that the underlying mechanism lies in the over-constraining of the model, and that the biaxial behaviour influences the identified Lankford coefficients. Overall, the three standard tests were successfully replaced with one deep-notched test performed at 45 • -User-defined virtual fields, selected based on the recommendations of [15], were used as an alternative to the SBVFs. It was found that the identified parameters were less accurate compared to the ones obtained with the SBVFs. The behaviour of the material was only matched along the orientation of the test used for characterisation (45 • ) and experimental forces in other directions could not be predicted correctly. -For Hill48 using total, as opposed to incremental, stress sensitivity fields to generate SBVFS, improved consistency of the identification. This is most likely due to incremental stress sensitivity fields filtering out the majority of the supplied data, especially once the entire ROI has yielded. On the contrary, when the total stress sensitivity fields are used, time steps with larger plastic deformation have more impact on the cost function and identification. -The number of load levels included in the identification did not significantly influence the values of identified material parameters for Hill48.
-Deep-notched tests were not rich enough to identify Yld2000-2D over all possible loading states. It was found that the identified and the reference yield surfaces matched well, but only in the regions represented in the test. When data from the bulge test was used to constrain the cost function, the identified model matched the reference very well over the entire stress space. -It was found that the commonly used Levenberg-Marquardt algorithm was not capable of finding the global minimum of the cost function. The sequential quadratic programming algorithm was used instead as it was robust enough to converge to the global minimum from a number of independent starting points.
The authors believe that the sensitivity-based virtual fields provide a significant step forward to calibrate non-linear models with the VFM, but there are still problems worth investigating. In particular, there are no guidelines on how much data should be used in identification (e.g. maximum strain, magnitude of strain increments). This should be studied rigorously, particularly employing the procedure of image deformation to account for the DIC parameters involved in the data processing [64]. Finally, the selection of the sample geometry still remains an open problem, and the results could certainly be improved further if the test was richer in terms of load paths. An important question to answer is how to design a better test that will contain sufficient information about the model and reduce the number of tests involved. Currently, numerical strategies are being developed, capable of optimising specimen geometries that promote heterogeneity and particular strain states [30,31]. They were applied to modify uniaxial tension specimens and although greatly improved heterogeneity of the test, struggled to produce data at biaxial tension thus had limited applicability to advanced constitutive models such as Yld2000-2D. Alternatively, a single modified cruciform specimen could be used, which produces data covering a large portion of the yield surface, including equibiaxial tension [63]. In future, the method presented in this contribution could be applied to those test configurations to identify models such as Yld2000-2D or Yld2004-18p. a Doctoral Training Grant studentship. Dr Davis also acknowledges support from the Leverhulme Early Career Fellowship. Dr Marek would also like to thank Dr Georges Limbert and Prof. Sandrine Thuillier for useful suggestions during his PhD examination. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Parameters F , G, H and N were identified with both of those methods and then used to calculate σ y 11 , σ y 22 , σ y 33 and σ y 12 for reporting using equation (7).

Yld2000-2D
The model is matched to 8 data points: r 0 , r 90 , r 45 , r b , σ 0 , σ 90 , σ 45 and σ b according to the procedure outlined in [6]. In essence, experimental flow ratios are fed to the model and the corresponding equivalent stresses and plastic flow ratios are calculated with equation (8) and the associated flow rule under appropriate stress states (uniaxial, biaxial). These observation are then compared against experimental measurements and the material parameters are adjusted such that the two match in the least-square sense.

Appendix B: Raw Data Report
This appendix contains detailed information about the data used in the report. Table 7 shows which images (available from the online repository) were used in the identifications carried out. Image offset is the number of the first image correlated in the DIC software; load steps taken indicate numbers of correlated images included in the identifications (notation used here is consistent with the Matlab notation).