Extension of the sensitivity-based virtual fields to large deformation anisotropic plasticity

The virtual fields method is an approach to inversely identify material parameters using full-field deformation data. In this work, we extend the sensitivity-based virtual fields to large deformation anisotropic plasticity. The method is firstly generalized to the finite deformation framework and then tested on numerical data obtained from a finite element model of a deep-notched specimen subjected to a tensile loading. We demonstrated the feasibility of the method for two anisotropic plasticity models: Hill48 and Yld2000-2D, and showed that all the parameters could be characterise from such a test. The sensitivity-based virtual fields performed better than the currently accepted standard approach of user-defined ones in terms of accuracy and robustness. The main advantage of the sensitivity-based virtual fields comes from the automation of virtual fields generation. The process can be applied to any geometry and any constitutive law.


Introduction
Accurate modelling of metal forming processes is of key interest to industries such as automotive. One of the main issues in simulating processes such as deep drawing of metal sheets is ensuring that the chosen constitutive model represents the material accurately. Many of the metallic materials used in this industry exhibit anisotropic properties due to texture induced during cold rolling which highly affect deformation of it during forming processes (e.g. earrings formation during deep drawing [20]).
One of the most popular anisotropic plasticity yield criteria is Hill48 [19], which in case of plane stress conditions requires four parameters, generally identified from three uniaxial tests performed in three directions: rolling (RD), transverse (TD) and 45 • to RD. In many cases, experimental results suggest that Hill48 is not capable of predicting biaxial yield behaviour accurately, and thus has limited applicability for forming predictions. Numerous models were proposed to capture the biaxial behaviour of sheet metals more accurately, such as: Yld89 [9], Stoughton's model of 2002 and further refinement of 2009 [43,44], BBC2000 [6], BBC2005 [5], Yld2000-2D [8] and Yld2004-18 [7]. Often, the usefulness of these complex models is limited by the significant effort required to accurately identify their parameters experimentally. In particular, many of these models require performing an additional biaxial test, such as bulge or equibiaxial tension on cruciform specimens, increasing the cost of the procedure. Therefore, there is a drive to improve testing techniques and a promising way to achieve this goal is to collect experimental data using more advanced methods, such as full-field measurements e.g. digital image correlation (DIC).
New tests can be designed in order to collect more data within a single run, compared to the standard methods. The use of full-field measurements makes it possible to choose complex geometries for the test specimens, introducing heterogeneous strain fields, thereby enabling the yield envelope to be probed at thousands of different stress states at once. One of the main challenges in such approach is to extract the material parameters from the collected data. Two of the most used inverse techniques capable of doing this are finite element model updating (FEMU) and the virtual fields method (VFM). It was demonstrated that these approaches can reduce the number of tests needed to fully characterize anisotropic models, in particular Hill48, or Yld2000-2D [11,13,18,21,35,38,39].
The virtual fields method is a very efficient technique for extraction of material parameters from full-field measurements. One of the main advantages of the VFM over the FEMU is that it is significantly faster in terms of the computational time. In fact, some authors reported that for their particular application the VFM was 125 times faster than FEMU [47]. This is especially important as the complexity of material models and the number of data points available from the measurements grow. Another advantage of the VFM is that it acts directly on collected data and no numerical simulations are required. As a result, the method can be integrated directly into a DIC platform making it more accessible to practising engineers. The method has already been applied to a range of materials and constitutive laws such as arteries [2], rubbers [17,45,46], composites [15], and metals [22,24,25,33,39].
One of the main challenges in the VFM is the choice of virtual fields. These are test functions that act upon the reconstructed stress fields to check for stress equilibrium. Their choice strongly affects the accuracy of the identification. Until recently, no structured method was available to generate high quality virtual fields for non-linear problems. Currently, the standard approach is to rely on user-defined virtual fields (UDVFs), using standard expansion bases such as polynomials or harmonic functions. The effectiveness of these user-defined virtual fields strongly depends on their choice, and requires the user to understand the method in depth to be able to select these fields in an informed way. Recently, a new approach for generating high quality virtual fields has been developed, leading to the so-called sensitivity-based virtual fields (SBVFs) [28]. They outperformed UDVFs in case of isotropic plasticity, and are generic enough to be applied to any constitutive model.
In this work, we have extended the SBVFs to the case of large deformation anisotropic plasticity, and demonstrated their feasibility to calibrate Hill48 and Yld2000-2D yield functions from a deep notched specimen subject to tensile loading.

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 motion of each material point can be described by a function x = φ(X, t), which maps the position of every particle in the reference configuration to the current deformed configuration. 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 [12].

The virtual fields method
Quasi-static equilibrium can be expressed in so-called 'weak form' in which it is enforced as a weighted average over the entire domain, expressed here in the current configuration B t , in absence of volume forces: where ∂B t is the boundary of B t , n is the outwards vector of ∂B t and σ is the Cauchy stress tensor. Equation 6, called the principle of virtual work (PVW), is satisfied for any continuous function u * (called virtual displacements) that is piecewise-differentiable. Both stress and test function (virtual displacements) are expressed in the current configuration in the case of Eq. 6. Fig. 1 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 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 The PVW can be alternatively formulated in the reference configuration. In that case, another stress tensor is defined, called the first Piola-Kirchhoff stress tensor P: Notably, the first Piola-Kirchhoff stress tensor is asymmetric due to the asymmetry of the deformation gradient. As a result, in the reference configuration, the PVW can be expressed as: where B 0 is the considered body in the reference configuration and ∂B 0 its boundary. This form is much more suitable for practical implementation in case of the proposed method, as the virtual fields U * , defined in the reference configuration B 0 , do not need updated virtual boundary conditions, as will become apparent later in the article. This approach has been used by most of the VFM community [36,38,39,41]. Noticeably, it was demonstrated that the current configuration formulation could be successfully applied to the case of hyperelasticity as well [2,23,32]. The VFM uses the PVW to identify material parameters from kinematic data and loading. Generally, kinematic fields are measured by means of full-field techniques such as DIC over the entire domain. This data is then used to reconstruct the stress field using a set of material parameters, denoted here as χ. The calculated stresses must be in equilibrium with the measured loading, which is enforced through either Eq. 6 or 8. As the correct values of the constitutive parameters are unknown at the start of the process, the stress field is first estimated with a guessed set of parameters and the equilibrium can be checked by means of Eq. 8: The material parameters are found through an iterative minimisation of the cost function, i.e. the correct material parameters produce a stress field that minimises the gap in the PVW. Since the full-field measurements provide spatially dense data, the integral of the stresses can be approximated by a discrete sum. Additionally, multiple load levels (time steps) have to be included to involve all the parameters of the constitutive model. Multiple independent virtual fields can be used to involve data in the cost function in various ways, which generally leads to better conditioning of the cost function, resulting in more robust minimisation and more accurate identification. Finally, a general form of the cost function can be expressed as: where W int is the virtual work of internal forces and W ext is the virtual work of external forces, S j is the surface area of j -th point and h is the thickness of specimen. Since in most of cases the measurements are only performed at the surface of specimen, some assumption about the through-thickness distribution of the mechanical fields must be made. Often, plane stress is assumed, provided the thickness of specimen is small in comparison to the other two dimensions and the loading is in plane. The measurements provide in-plane kinematic quantities and the out-of-plane stresses are considered negligible (σ 13 = σ 23 = σ 33 = 0). However, when the PVW is formulated in the reference configuration (8), the 2D Cauchy stress must be pulled-back from B t to B 0 , as shown in Eq. 7. In reality, the gradient of deformation is fully 3D and so some additional assumptions must be made. Assuming that for a thin specimen the out-of-plane shearing is negligible (F 13 = F 23 = F 31 = F 32 = 0), the Jacobian (det(F)) can be now expressed as: The in-plane values are measured, however, the out-of-plane term is still unknown. It can be approximated by calculating the out-of-plane strain during stress reconstruction using the elasto-plastic constitutive law: Here, the Hooke's law was assumed for the elastic part and the isochoric flow for the plastic part of the strain increments. This strain can be then used to calculate the missing component: It should be noted that if the isochoric flow assumption becomes questionable, back-to-back camera systems can be used to determine an average value of ε 33 , as shown in [14]. Finally, with known Jacobian it is possible to pull-back the Cauchy stress to the reference configuration.

Sensitivity-based virtual fields
The test functions (virtual fields) in Eq. 10 are arbitrary and must be selected before the identification is conducted. The selection of virtual fields has a significant impact on the accuracy of the identification. These functions influence the amount of error introduced to the cost function by selecting which data points, and with what weight, are introduced to the cost function. The main difference between virtual fields arises in the way they propagate experimental noise. In linear elasticity, an automated procedure has been published in 2004, relying on the minimization of the impact of noise on the identified parameters, i.e., finding the virtual fields leading to the maximum likelihood solution for a given basis of functions to expand the virtual fields [3]. This is now routinely used by the VFM community and also implemented on the commercial DIC/VFM platform MatchID [29]. An attempt at extending this to non-linear laws, namely, isotropic elasto-plasticity, was published in 2010 [34]. The idea there was to use a piecewise linear definition of the virtual fields based on the tangent stiffness matrix. The method did improve results but was found to lack flexibility as it required an expression for the tangent matrix. Also, for non-linear models where strains are generally larger, sensitivity to noise is not necessarily the most relevant criterion to select virtual field.
Recently, a new type of virtual field for the non-linear laws was proposed [28]. They are automatically generated during the identification procedure with very limited user input. These fields, called sensitivity-based virtual fields, are based on the reconstructed stress field and so, easily and automatically adapt to any geometry and material model. They were shown to outperform the user-defined virtual fields for isotropic small-strain plasticity and seemed very promising for more complex problems. They were also shown to outperform the tangent-matrix fields from [34], though only marginally for this particular case.
The main idea behind the sensitivity-based virtual fields is to find areas during the test, both in space and time, where the information about each parameter is contained. A separate virtual field is constructed automatically for each constitutive parameter which allows the cost function to represent each parameter with maximised sensitivity. In order to locate areas where the information is encoded for the i-th parameter, an incremental stress sensitivity is calculated. By perturbing the value of the i-th parameter during the stress reconstruction for a given set of current parameters, a change in stress field is noted, highlighting areas where the parameter is active in influencing the stress. This information is encoded in a map called stress sensitivity: where δχ i is a perturbation of the i-th parameter, typically −0.2χ i ≤ δχ i ≤ −0.1χ i . The negative sign is taken to expand the VFs over points currently not active, as opposed to penalizing points that just became active (e.g. in the case of yield stresses). Furthermore, incremental stress sensitivity maps are found as: These sensitivity maps, either incremental (15) or total (14) can be used as virtual strains to provide relevant weight to the stresses in the PVW equation. This incremental approach was proposed in order to effectively decouple the influence of yielding-related parameters from hardening parameters which act at different time scales. The incremental sensitivity maps cannot however be directly used in the PVW as the corresponding virtual displacements, needed in the PVW, are unknown. However, it is possible to construct virtual displacements such that the resulting virtual strain fields 'look like' the incremental sensitivity maps. This can be achieved by performing a least-square match, under some constraints, between virtual strain fields and incremental stress sensitivity.
An effective way to solve the matching problem is to construct a virtual mesh, which consists of virtual nodes connected through virtual elements. The virtual mesh defines the virtual displacements at the nodes and interpolate them within an element using classical FE shape functions. The nodes are connected with 4-node quadrilateral elements and the interpolation is done using standard bi-linear shape functions (N). The virtual fields can be calculated using the values of virtual displacements at each of the node and the spatial derivatives of the shape functions. For every point of measurement, the virtual strain fields can be found as a function of virtual displacements of the neighbouring nodes. This function is linear with respect to the displacements and can be expressed using the strain-displacement matrix B el , which relates virtual displacements at the nodes to the virtual strains at the considered data point: where U * is a vector containing the virtual displacements of the element nodes and ∂U * ∂X is a vector containing values of virtual strains 1 at the data point. It is chosen here to express the virtual fields in the reference configuration, as a results B el does not change with deformation. It can be expressed as: For every point, matrix B el can be generated and then assembled into a global strain-displacement matrix (B glob ) which relates virtual displacements from all virtual nodes to the virtual strains at all data points. The global straindisplacement matrix has then to be modified to account for virtual boundary conditions. Often, it is necessary to constrain the virtual displacement in the direction of applied loading to be constant across the loading edge, so that only the resultant load appears in the PVW equation and not its (unknown) distribution. This simplifies the traction contribution to the PVW (6,8) to become: (18) where F load is the total force applied. Imposing these constraints into B glob , a modified global strain-displacement matrix is found,B glob . In order to solve the least-square problem, a pseudo-inverse of this matrix is found which can be used to generate the corresponding virtual displacements: These virtual displacements now produce virtual strain fields that 'look like' the incremental stress sensitivity maps and obey the necessary virtual boundary conditions. The virtual strain fields are finally found with the following formula: Note that the construction of the sensitivity-based virtual fields must be performed at every time step. However, as mentioned before, if the reference configuration is chosen for the PVW, the matrix B glob is assembled only once for the entire identification.
Computing stress sensitivities significantly increases the identification time, as it virtually doubles the number of necessary stress reconstructions. To improve the computational efficiency, a selective updating scheme can be employed. Recall that any continuous virtual displacement fields constitute a valid choice, including the sensitivity-based virtual fields based on incorrect (e.g. initial) parameters. Effectively, these can be used to put the minimisation algorithm in the neighbourhood of the solution without updating them, but carrying across the iterations. As the algorithm converges, the virtual fields can be updated with parameters much closer to the correct values, saving many stress reconstructions and significantly improving computational efficiency.
Finally, in order to balance the contributions from each virtual field, a scaling is introduced. The virtual displacements are scaled by a factor dependent on the current internal virtual work contributions (W int ). For each iteration, W int is calculated (10), and then sorted according to the absolute values over all time steps. The scaling parameter is calculated as a mean out of the top x th percentile of the sorted values. This ensures that virtual fields contributions associated with each parameter are of similar orders of magnitude.

Numerical simulations
The standard way to test anisotropic materials is to conduct tensile tests on dog-bone specimens cut at different angles to the rolling direction of the sheet (typically 0 • /45 • /90 • ). If the material model under inspection includes a parameter related to biaxial yield stress, an additional test is required, either a bulge test or equibiaxial tension on a cruciform specimen. The major limitation of this approach is that a single test provides only one data point on the yield locus and many tests are needed to match the yield surface.
An alternative is to run a test with enough stress heterogeneity to identify all necessary parameters at once. Some of the heterogeneity in the tensile test can be obtained by means of material orientation, geometric features, and loading. Rossi et al. [39] proposed a test on a deep-notched specimen under tensile loading capable of identifying the Hill48 model using a single specimen. The test is replicated here and combined with the sensitivity-based virtual fields to test their applicability to large strain anisotropic plasticity. This does not mean that this test is optimal in any way, but it can serve as a clear comparison on how VFs selection impacts the identification. Many different geometries have been proposed in the literature to produce heterogeneous states of stress and strain, it is beyond the scope of the present paper to investigate this. However, future work will look at specimen optimization, in the same spirit as for composites testing in elasticity [16].

FE model
The test proposed by Rossi et al. involved tensile loading applied to a flat coupon with two circular notches, which introduce heterogeneous deformation [39]. It was simulated in Abaqus (v. 6.13) to generate synthetic data which were then used to test the identification algorithm. The geometry of the specimen is presented in Fig. 2. The mesh density was chosen according to a convergence study. The thickness of the plane stress elements was chosen as 0.74 mm, similar to that used in [38] and typical of thin anisotropic metal sheets for the automotive industry. The bottom edge was fixed, while a constant vertical displacement of 6.75 mm was applied across the top edge, which was additionally fixed in the lateral direction to simulate the effect of the grip of a test machine. The initial material orientation was defined by specifying the angle (θ ) between the rolling direction and the horizontal axis of the model, indicating the principal material axes (Ξ, H ) needed to describe the anisotropic properties (see Fig. 2). Overall, total vertical strains of about 40% were obtained with the model.

Constitutive models
In the VFM, the stress field is reconstructed explicitly from the kinematic measurements through an assumed constitutive law. In this work, two different large strain plasticity models were considered: Hill48 and Yld2000-2D.
The elastic response was modelled with Hooke's law extended to finite deformation: where Δσ is the rate of change of the Cauchy stress, D is the elastic operator for plane stress and Δε e L is an increment of the elastic part of the Hencky strain. Here, we assume that both Young's modulus E, and Poisson's ratio ν needed to construct D are known. They can be identified from full-field data using the noise-optimised virtual fields in elasticity [4], considering only the initial elastic steps.
In order to ensure the objectivity of the model, a corotational frame was adopted. The reference frame at each material point rotates with the material. Fig. 1 presents the corotational and material frames in both reference and current configurations. Initially, the corotational frame (i, j ) is co-linear with the global frame. The angle between the material coordinate system (Ξ, H ) and the global reference frame is known. This angle can be used to construct a tensor (R mat ) rotating the local coordinate system to the material one. Due to deformation, the corotational frame rotates with R and in the current configuration is denoted (1,2). The angle between the current material coordinate system (ξ, η) and the corotational frame is fixed and related to R mat . The strain increment in the current material frame (ξ, η) can be expressed in terms of the strain increment in the global frame (i, j ): More complex strategies could be employed to account for the evolution of material texture such as that presented in [31]. They would have to be implemented in the constitutive routine used to reconstruct stress (Fig. 4). Consequently, an inverse rotation could be performed in order to represent the reconstructed stress tensor in the global coordinate system, in which the PVW is expressed. The strain increments in the global coordinate system are simply calculated as the difference between two consecutive total strains: For each of the two plastic models investigated, Hill48 and Yld2000-2D, an associated flow rule was assumed, as well as an additive decomposition of strain increments: where Δε e L is the elastic strain increment and Δε p L is the plastic strain increment.
A yield criterion, in general, can be written as: where σ eq is an equivalent stress formulated differently for each plasticity model and σ y is a yield stress evolving according to a hardening law. A shared feature of both models is that they are formulated in terms of Cauchy stress expressed in the current material frame (ξ, η in Fig. 1).
Hill48 is a popular model for anisotropic plasticity [19]. When plane stress is assumed the equivalent stress can be expressed as: 2 12 . (26) This model depends on 4 independent parameters defining the anisotropy. A suitable way to express the governing parameters is to relate them to mechanical quantities measured in an experiment. In this paper the plastic potentials (R ij ) are used to define the criterion: where R ij = σ y ij σ 0 . Note that σ y 11 and σ y 22 are the yield stresses identified in planar uniaxial tests conducted at 0 • and 90 • respectively and σ y 33 is the through-thickness yield stress. Finally, σ y 12 is the yield stress identified under pure shear. Furthermore, it was assumed that the yield stress in the hardening law is equal to σ 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 model.
Although the model is popular in literature, it suffers from poor performance when used in context of sheet materials subject to biaxial loading [30,39,44]. This is mostly due to the quadratic nature of the law which does not represent real materials accurately. It was found that nonquadratic models such as Yld2000-2D or Stoughton2009 predict the behaviour of sheet metals such as steel or aluminium more accurately [44].
Yld2000-2D [8] was developed strictly for plane stress conditions 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 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 , and can accurately represent the behaviour both in simple tension as well as biaxial loading.
Two different hardening laws were adopted here: linear (32) and the power law hardening (33): Linear hardening is defined with two parameters: initial yield stress σ 0 and hardening modulus H . The nonlinear power law includes 3 parameters: initial yield stress σ 0 , hardening modulus H , and an exponent n. The equivalent plastic strainε p is integrated over the history of deformation by means of summing the equivalent plastic strain increments obtained at each increment: The reference parameters used to generate the models are presented in Tables 1 and 2. For each of the constitutive models a routine was produced integrating the constitutive equations using an implicit scheme with a radial-return algorithm, returning the state of stress and out-of-plane strain at each time step, with both being rotated back to the global coordinate system. For Hill48, 6 data sets were generated in total, considering both hardening models at the three different material orientations (30 • , 45 • , 60 • ). Due to its computational intensity, only one model was considered for Yld2000-2D, simulating linear hardening with a material orientation of 45 • . While in this work the constitutive modelling was kept simple, the proposed methodology is valid for any chosen material model. In practice, the investigator supplies a constitutive law to be identified (see Fig. 4). The model can be arbitrarily complex, as long as it is capable of reconstructing the stress field based on the measured kinematic data (e.g. deformation gradient) and internal state variables that can be carried over between different load levels. By choosing more complex models that account for e.g. multiplicative decomposition of deformation gradient [27], hyperelasticity [42], anisotropic hardening [10] or even complete anisotropic elastoplasticity [40] more complex material description could be reached, but this has not been considered in this work.

Data processing
Following the FE simulations, raw displacements and positions of data points were extracted from Abaqus and exported into Matlab (2016a). The cloud of points was then interpolated onto a rectangular grid of 409 × 349 data points, corresponding to the data density typical of a DIC measurement. The region of interest was then trimmed to a grid of 362 × 202 points spanning the region of interest (ROI) (Fig. 2), which produced 57,392 data points. In the case of Yld2000-2D model, to reduce the computational time, a coarser data grid was used: 214 × 119 producing 20,024 points in total.
The displacements were then corrupted with a Gaussian white noise, with standard deviation of 0.3 μm, representative of a real experiment. For each data point, displacements were smoothed using spatial and temporal filters, as typically done for experimental data. Temporal smoothing was performed using the Savitzky-Golay method with a polynomial order m temp and a window size of w temp . This was complemented with a spatial smoothing using a Gaussian filter of standard deviation σ spat , with the window size adjusted to be the maximum odd number smaller than 3 × σ spat × 2. While temporal smoothing was performed on all 400 time steps obtained from Abaqus, only some of the temporal data points were passed to the identification procedure. The main reason for that was to decrease the computational effort, but also this increased the size of the strain increments mitigating the impact of strain noise on the error in the stress reconstruction. The time steps used for the identification are graphically presented on the global force-displacement curves in Fig. 3.
For each model, a different number of time steps was taken. In the case of Hill48 with linear hardening, the calculations were relatively fast, thus having many temporal points was not a problem. Ultimately, 83 frames were considered. In the case of the power law hardening, during deformation a plastic instability occurred and the strains began to localize in a very small band, causing a geometrical softening. This data was discarded, as shown in the figure.
Overall, only 60 frames were used, with a maximum plastic strain of about 30%. Finally, for Yld2000-2D only 27 time steps were used, as the constitutive law itself is much more computationally demanding, and having more points would just lead to inconvenient computational times, unnecessary at this first validation stage.
After smoothing and down-sampling were performed, the deformation gradient F was calculated using central finite difference, which was then used to calculate V with Eq. 4, R with Eq. 3 and ε L with Eq. 5. The set of the three quantities (F, R, ε L ) for all data points constitutes the kinematic data used to identify the material parameters.
The kinematic data was then used as an input to an in-house Matlab code implementation of the VFM. The stresses were calculated using the deformation data and an initial guess for the material parameters. The SBVFs were calculated during the stress reconstruction process which allowed the values of the cost function to be calculated with Eq. 10. The material parameters were refined iteratively by means of the Matlab function fmincon minimising (10) with a sequential quadratic programming algorithm (SQP). Starting points were selected randomly from the interval between 50% and 200% of the reference values using a random number generator. Two starting points were used for   each data set in order to ensure that the identified parameters corresponded the the global minimum. The identification algorithm is summarized in Fig. 4 in the form of a flowchart. Finally, the identified parameters were compared against the reference values to quantify the accuracy of procedure. The SBVFs were updated when the first-order optimality [1] fell below certain threshold (1 × 10 −4 ), indicating convergence of the procedure. The virtual fields were recomputed and stored in memory, and the threshold was scaled down by a factor of 2.3 to allow further refinement. Sometimes, after selective updating was performed, the value of cost function would increase, creating an apparent local minimum terminating the minimisation. To prevent this, the value of cost function was offset below previous iteration when the update was performed. This scheme leads to an increasing rate of updates as the solution converges, providing valid values of virtual fields at the optimum. In total, about 10-15 virtual fields were computed throughout the identification consisting of more than 100 iterations, with the number controlled indirectly by the two parameters (first-order optimality threshold and threshold refinement parameter).

User-defined virtual fields
This works aims at extending the SBVFs to large deformation anisotropic plasticity. As mentioned earlier, the specimen design was previously proposed by Rossi et al. [39] who employed UDVFs to identify Hill48 parameters using a single test, as well as Yld2000-2D parameters using a combination of three tests. In this work, their user-defined virtual fields were used as a benchmark to test the suitability of the SBVFs. They showed that suitable VFs for the test consisted of a combination of three virtual fields, defined with the following virtual displacements: The virtual displacements are defined in the coordinate system presented in Fig. 5. The virtual fields are constructed in a way to include all stress components in the cost function.
It is worth noticing, that only the first virtual field (35) includes the virtual work of external forces. This impacts the balance of contributions from each virtual field to the cost function. Traditionally, the residual is scaled by the virtual work of external forces to provide a dimensionless value Fig. 4 Flowchart representing the algorithm for material identification. There are three major parts: generation of raw data, processing of raw data and material identification with the VFM [33]. This was not possible here, as two of the virtual fields do not include the contribution of the external forces. To overcome this problem, scaling by the maximum value of the internal virtual work was introduced, which normalises the peaks of the internal virtual work to 1. To balance all three fields even further, the residual coming from the first virtual field was scaled by a factor 500, which was found to produce the best results. This workaround scaling shows that the manually defined virtual fields lack generality, and must be individually tailored for each application, which is a significant disadvantage.

Error quantification
As the number of material parameters in the constitutive models increases, quantifying the accuracy of identification becomes a challenging task. Often, these models are defined with parameters that lack physical meaning, and on its own have limited impact on the model outcome, which is driven by the compound action of all of the parameters, as for Yld2000-2D. It is therefore important to establish tools to compare different sets of identified parameters in order to quantify the accuracy of the identification meaningfully.
Because of the rather large number of parameters and their individual lack of physical meaning, a direct comparison between reference and identified parameters on a one to one basis is not always relevant to draw meaningful conclusions. The apparent uniaxial yield stresŝ σ as a function of material orientation was found to provide a convenient and physical-based way of evaluating identification errors. Assuming a yield criterion of the form of Eq. 25 and a uniaxial state of stress at angle θ with respect to the material coordinate system, the stress state in the material coordinate system can be expressed as: This stress can be introduced into the expression of σ eq and the criterion can now be algebraically transformed into: Using Eq. 39, it is possible to reconstruct the variation of the yield stress with angle θ and equivalent plastic strain ε p . Such map can be calculated for both the reference and identified parameters, so that they can be meaningfully compared. The error can be quantified by looking at the map of relative error: The map can be used to investigate whether the error is associated with the variation of yield stress with material orientation and/or hardening. Alternatively, the mean of the error map,r, gives a single numerical value for the accuracy of the identification, later referred to as the global error: This approach enables the reference and identified sets of parameters to be compared in a more meaningful way as opposed to parameter by parameter and furthermore, it also allows for meaningful comparisons across different material models.

Validation of the method
As a first validation of the methodology, the six raw FE data sets generated for the Hill48 model were used. Since the data did not contain any noise, no smoothing was used.
The following parameters were selected for generation of the sensitivity-based virtual fields: δχ i = −0.15χ i , virtual mesh of 14 × 14 elements, and the scaling parameters based on the mean of top 30 th percentile of the internal virtual work terms. As shown before, the choice of the parameters has a limited impact on the identification and mostly affects the random portion of the error [28]. The results of the identifications are presented in Table 3. Clearly, for linear hardening, all tests were successful at identifying the material parameters. Although marginal here, it is worth noting that the mean error for sensitivity-based VFs is already consistently smaller compared to the user-defined ones. The single tests performed at either 30 • or 60 • were only performed usig exact data and the linear hardening. As they are unlikely to be practical when experimental and modelling errors are introduced, they were disregarded for the more complex cases. For the power law hardening, the identified values were also very close to the reference. With the exception of the combined 30 • +60 • tests using the userdefined VFs, the mean error was consistently below 2%, showing that the methodology adopted here has a potential of identifying Hill48 parameters using a single test. The combined test was adopted in order to obtain more data over uniaxial states of stress, in comparison to the 45 • test. The latter, probes the yield envelope under combined shear and normal loading, with very little data points containing information about pure σ 11 or σ 22 behaviour. By including both 30 • and 60 • tests, more information is available about these regions, as well about shearing behaviour enabling identification all of the parameters.
Scaling of UDVFs is of crucial importance. Without introducing the scaling discussed in Section "User-defined virtual fields", the UDVFs lead to large identification errors. In particular, the yield stress is correctly identified only near the material orientation angle, while the remaining  (Fig. 6a and b). The scaling significantly improves the identification: see the difference between Fig. 6b and c orientations are characterised incorrectly, as shown in Fig. 6a and b. The reason for this is most likely due to the first set of UDVFs (35) being overrepresented in the cost function. As a result, the vertical component of the stress field represented in the global frame is the major source of information, leading to insensitivity of the cost function to the other two stress components. The identification can be improved by introducing the scaling which restores the balance between the three virtual fields. In that case the anisotropy is recognized much better, as shown in the difference between Fig. 6b and c.

Effect of noise in the data
Having validated the methodology, the effect of noise on the identified parameters was studied. Since in practice, spatial and temporal smoothing are always implemented, we also implemented this here to simulate the experimental process more realistically. Without any smoothing, the difference between UDVFs and SBVFs was even larger, but this was thought to be a somewhat unfair comparison. First, a sensitivity study was conducted to determine the optimal smoothing parameters. The study was performed using fixed settings, i.e. sensitivity-based virtual fields, linear hardening and the same level of noise in each run. In total, 30 copies of noise were processed. The spatial smoothing parameters were chosen under fixed temporal smoothing, likewise the temporal parameters were chosen under fixed spatial smoothing. As presented in Table 4, the smoothing has a limited effect on the magnitude of both systematic and random errors (global error value and its spread respectively), which suggests that the choice of parameters is not critical for the identification here. Generally, the stronger the smoothing, the smaller the random portion of error, at a cost of increased systematic error. However, it was found previously that significant noise can cause spurious elastic unloadings, which strongly affect both errors [26]. As a result, some smoothing settings can improve both errors at the same time. In order to balance systematic and random errors, spatial smoothing with a window of 9 was selected, combined with 11 points temporal smoothing using 3 rd order polynomial. The selected smoothing parameters were used to study the effect of the virtual fields on the systematic and random errors. The identification was run on 30 copies of noise, as this gives enough statistical representation to establish the random part of the error. The results are presented in Tables 5 and 6 for the linear and the power law respectively. The parameters are expressed relative to the reference values and the uncertainty represents a coefficient of variation expressed in % The parameters are expressed relative to the reference values and the uncertainty represents a coefficient of variation expressed in %

Linear hardening
In the case of linear hardening, the sensitivity-based virtual fields provide more accurate identification, however with about three times larger random error compared to the UDVFs. Although individual parameters exhibit errors as large as 5%, the global error is significantly smaller, indicating compensation between parameters. Remarkably, the overall error is very small, for the level of white noise added to the displacements. It must be noted that when no smoothing was used on the data, the accuracy was poor and the global error was about 20%. By replacing a single test (45 • ) with two combined tests (30 • + 60 • ), the random part of the error is reduced by a factor of 4. This was expected as more information is used, which averages out noise more effectively. The effect on the systematic error is not as notable; for the sensitivity-based VFs the mean error is slightly reduced, while for the user-defined VFs it increased a bit, as shown in Fig. 7a.

Power law hardening
The results are markedly different for the power law hardening. Firstly, the effect of noise is significantly larger compared to the linear hardening case, as shown in Fig. 7 (note the difference in scaled between Fig. 7a and b). The mean errors for the 45 • tests are between 4 and 5%, compared to about 2% for the simpler hardening law (Fig. 7). While level of error is still satisfactory, it must be noted that it represents a lower bound as the simulation of experimental uncertainties remains very basic. Unlike the linear hardening case, the random error was the smallest when the SBVFs were used and the smallest systematic error was obtained with the UDVFs. Significant improvement was found when the single test was replaced by the combined tests, both in terms of systematic and random errors. This indicates that the additional test contributes significant data to the cost function. Since the smoothing parameters were chosen for the linear hardening model, it is possible that a different combination of parameters would provide smaller errors, and differentiate between the two virtual fields types in a different way. While the single mean error term makes it straightforward to compare the overall accuracy of the identification the influence of errors on the yielding and hardening parameters can be more readily understood by examining the error maps, proposed in Section "Error quantification". In Fig. 8 the SBVFs provided the larger errors at 45 • , even though the single test at 45 • contains direct information about yielding at this orientation. For the UDVFs the error was the lowest for angles close to the specimen orientation, for the single test at 45 • and for the combined tests at both ends of the angle spectrum, as expected from a mechanical perspective. A possible explanation is the overconstraining of the sensitivity-based virtual fields. As four independent stress sensitivity components are mapped onto only two virtual displacements maps, the resulting virtual fields are not reproduced perfectly. This appears to be especially true for the shearing components as they are cross-derivatives of the calculated virtual displacements. This might enhance the error on the shearing yield stress, observed as a slightly larger error at 45 • . Further study needs to be done to confirm  [39], which was time-consuming. The SBVFs provide nearly equivalent results using a systematic procedure without any trial and error, resulting in a much faster and more rigorous process.
The quality of the virtual fields can be indirectly quantified by means of the curvature of the cost function at the minimum. Theoretically, the more curved the cost function, the more stable the identification. Additionally, where H V F is the Hessian matrix for a given virtual field, Φ ident is the value of cost function at minimum, and ∂ 2 Φ V F ∂χ i ∂χ j is the second derivative of the cost function with respect to the material parameters. The balance of the terms is mostly dominated by the mechanics of the test, however it can be improved by a good choice of virtual fields. In order to investigate how the balance is affected by this choice, it is convenient to look at the Hessian scaled by the largest term. The curvatures of the cost function for the linear hardening case are presented in the figure below: By comparing the two Hessians ( Fig. 9) it is apparent that the SBVFs do slightly better job at balancing σ

Yld2000-2D identification
In the case of the Yld2000-2D model, the main goal was to explore the comparative robustness of the UD and SB VFs when the model is richer in parameters to be identified. Only exact data from the 45 • test were used, as it is already a good example of the UDVFs underperforming. The identified parameters are presented in Table 7. It is worth noting, that in the case of Yld2000-2D, the anisotropic parameters do not have any obvious physical meaning, thus comparing them on a parameter to parameter basis is not so relevant to draw conclusions as discussed in Section "Error quantification". The distribution of initial and final yield stresses are presented in Fig. 10 and the identified yield loci are shown in Fig. 11. The superiority of the SBVFs is spectacular. The distribution of anisotropy is very well identified and the main source of error comes from the hardening modulus, which is consistent with the trend observed for Hill48. In the case of the UDVFs, the identification was highly inaccurate. The parameters found were significantly different from the reference, and four out of eight reached identification constraints (set to 50% and 200% of the reference). Most likely, this is due to the generic nature of the fields. Yld2000-2D is a more complex model  Fig. 12.
In both cases, there is a strong sensitivity to the shearing components (α 7 , α 8 ) and a moderate one for α 4 . When the UDVFs were used, there was very little sensitivity to the remaining parameters, except for the cross-correlation between α 7 and three other parameters: α 3 , α 5 and α 8 . In contrast, for the SBVFs, many more parameters were active. There was a cross-correlations between α 7 and all other parameters, as well as small to moderate sensitivities for all Reference SBVF UDVF Fig. 11 Initial shape of the yield stress surface for the Yld2000-2D model. The identified curves were generated using parameters identified in noiseless 45 • test data. The outline is drawn at 0 shear stress α parameters, showing that most of the material parameters were active during the test and the SBVFs were capable of balancing their contribution.

Computational efficiency
One of the important practical aspects of the type of methodology presented here concerns computational efficiency. Information on computing times are often left out in publications on this topic and therefore, it is difficult to compare computational efficiency with competing techniques motivating our decision to report this information.
The identification was run on a standard PC with an Intel Core i5 processor (3.20GHz) and 24 GB of RAM memory. The stress reconstruction routines were coded in Matlab, however they were automatically translated to C language and called as mex file using Matlab coder tool. This procedure improved the efficiency of the stress reconstruction routine, which is the most computationally demanding process. The approximate running times of identification for Hill48 model can be found in Fig. 13. In case of Yld2000-2D the time needed to obtain the results were 122 and 278 hours respectively for the UDVFs and SBVFs.
In the case of Hill48 model, the identification took approximately six hours for linear hardening and between eight and ten hours in the case of power law hardening, due to the additional unknown parameter. The SBVFs are not significantly slower compared to the UDVFs, and the difference is mostly due to the additional stress reconstructions needed to calculate the incremental stress sensitivities. The number of these reconstructions can be controlled by the number of times the SBVFs are updated, and so the difference could be brought down even more. It should also be emphasized that these computing times would be significantly reduced using a compiled programming language instead of Matlab.
In the case of the Yld2000-2D model the identification times were much longer, mostly due to very slow stress reconstruction procedure, which on average took ten minutes for the chosen data density. Additionally, as 10 unknown parameters were sought, the calculations of the gradients in the minimisation problem call for many stress reconstructions, further increasing the disparity between Hill48 and Yld2000-2D. As a comparison, the FEM model used to generate the data, took 4 hours to complete.
Although the running times were very high, there are many ways in which they could be improved. First, and as stated before, efficient implementation in a fast language would lead to a significant improvement. Second, replacing the implicit stress reconstruction algorithm with a direct method such as the one used in [39] would make the stress reconstruction faster, especially in the case of Yld2000-2D for which computing derivatives for the implicit scheme is a very computationally intensive process. The direct method is valid only for large deformation data (with plastic flow well established), however it could be coupled with the implicit algorithm for elasto-plastic transition. It was found that when the SQP algorithm was replaced with the Levenberg-Marquardt algorithm, the computations were up to 10 times faster, indicating that choosing a proper tool for minimising the cost function is crucial. For instance, the identification of the Yld2000-2D model using the SBVFs was reduced from 278 hours with SQP to merely 26 with Levenberg-Marquardt. The major advantage of the Levenberg-Marquardt algorithm is that it requires significantly fewer iteration to find a minimum, compared to the SQP algorithm. In the case of Hill48 identifications the former converged in approximately 8 iterations, compared to 130 of the latter. However, the efficiency of the Levenberg-Marquardt algorithm heavily relies on the initial guess which must be close to the solution to obtain fast convergence. This is not required by SQP which can find the solution from any starting point.
Choosing an optimal data density remains an open problem. In this work, about 60,000 spatial data points were used for the identification problem. However, this could potentially be reduced. For DIC measurements, the number of data points could be effectively controlled by means of the stepsize. Additionally, the number of load steps taken into consideration could be varied. It is worth noting that if temporal smoothing is performed on the entire collected data, even the time steps not used explicitly in the identification affect the outcome, as the information is passed through the temporal filter. The limiting factor for the temporal resolution is related to the stress reconstruction algorithm, i.e. in the case of implicit algorithms the larger the strain increments, the larger the reconstruction error. More studies are needed on the optimal number of data points for accurate reconstruction of material parameters. To define optimal spatial and temporal sampling leading to acceptable systematic and random errors, the impact of DIC parameters, and smoothing would need to be assessed with a synthetic image deformation procedure as in [37].

Conclusions and future work
In this work we extended the sensitivity-based virtual fields, originally proposed in [28] to large deformation anisotropic plasticity, which is the main novelty of this contribution. We tested the performance of the fields using a deep-notched tensile test, already used in that context by Rossi et al. [39], and found that the proposed virtual fields can be successfully applied to the problem.
It was found that for the Hill48 model both UDVFs and SBVFs were capable of correctly identifying the parameters from a single test. It must be noted that the UDVFs were developed using a trial-and-error approach by the lead author of [39]. In that context, the systematic procedure of SBVFs seems to be especially appealing as it removes the need of an informed input by the user to arrive at correct parameters. This is especially important with even more complex models such as Yld2000-2D where deriving appropriate UDVFs is a challenging task.
The main advantage of the new fields comes from the automatic generation procedure, with limited input from the investigator, involved only in setting virtual mesh density and scaling parameters. It was found before that these parameters have very limited effect on the overall identification and can be chosen a priori without a significant impact on the parameter values [28]. As a result, high quality virtual fields are generated for any material model, regardless of the test configuration. This opens up possible implementation in a user-friendly VFM software for non-linear model identification, like the MatchID DIC/VFM platform [29]. We demonstrated the effectiveness of the method on both Hill48 and Yld2000-2D, with the latter especially successful in comparison to the standard approach.
The identifiability of Hill48 model from a single test has been already established before [13,39]. Interestingly, the results suggest that a single test performed at 45 • may contain enough information to characterise the Yld2000-2D criterion as well. In fact, if confirmed experimentally, it would give an exciting alternative to the standard test protocol involving three uniaxial tests and one biaxial, significantly reducing the experimental effort to characterise a material. This would be possible due to the ability of the SBVFs to identify in space and time when each parameter is active and focus exclusively on those regions when identifying the parameters.
The method is currently being validated experimentally on an automotive DC04 steel alloy. The results from both UDVFs and SBVFs will be compared to the parameters identified with the standard multi-test protocol to confirm whether the Yld2000 criterion can indeed be identified from a single heterogeneous test, which would be a significant step forward to reduce identification costs and time scales.
Although a relatively simple material model was employed in this work, the presented method is general and can be used with any constitutive model. There are no limitations on the complexity of material model used within the VFM framework, given that it reconstructs stress field from measured kinematic fields (deformation gradient) and some internal state variables (that can be resolved by considering history of deformation). It must be stressed that usually the more complex the model, the more material parameters must be identified experimentally. The method proposed here allows for complete identification of material parameters given the experiment contains sufficient information. To the authors' best knowledge there is no systematic way of assessing the level of information that a test contains given a constitutive model. Investigating this in future would be certainly of importance for designing better experiments.
The new route to virtual field selection demonstrated here has potential in many applications, in particular for non-linear models with large numbers of parameters. An obvious extension would be hyper-visco-elastic models. The VFM has been applied to such materials in the past, see [17] for instance, but always with UDVFs, limiting the complexity of the considered models. Another area of interest concerns transient dynamic tests to identify the high strain rate elasto-plastic response of materials. Finally, now that a systematic route has been clearly identified to generated virtual fields automatically for non-linear laws, the problem of test optimization can be addressed. This has been studied for linear elasticity in the past, thanks to the availability of the noise-optimized virtual fields from [3], and can now be addressed for non-linear models using a procedure similar to that in [16].