Enhanced features of a constitutive equation gap identification method for heterogeneous elastoplastic behaviours

To identify mechanical properties in heterogeneous materials, the local stress fields have to be estimated. The recent developments in imaging techniques allow reaching precise and spatially dense kinematic fields (e.g. displacement, strain ...). In this paper, an iterative procedure is used to identify the distribution of elastoplastic material parameters and the local stress fields. The formulation and the principle of the method are briefly presented while attention is paid to check its reliability and efficiency on finite element simulation data as reference full-field measurements. The method is also applied to noisy measured displacement fields to assess its robustness.

In this work, we extend and adapt the approach developed in [24] to identify the constitutive laws and their mechanical parameters for heterogeneous materials. Since the method proposed in [24] is based on Airy functions, the approach is limited to simple geometries and regular meshes. This limitation is here removed and any geometry can be addressed. Moreover, the initial work [24] was limited to elastoplastic behaviours with a linear hardening. For sake of simplicity, the present paper also focuses on linear hardening but it is now straightforward to deal with any kinematic hardening law. The last improvement presented here concerns the identification of the yield stress and of the hardening modulus. The formulation was modified in order to allow the simultaneous identification of the two plastic parameters: the initial requirement of a plastic zone whose size remains constant on two successive load steps is no more needed.
The new formulation thus allows to deal with multilinear hardening behaviours on complex geometries and several load steps. The simultaneous use of several load steps for the identification significantly decreases the sensitivity to measurement noise.
The class of models that we have in mind belongs to J2 elastoplasticity with hardening. The Constitutive Equation Gap Method (CEGM) originally used as an error estimator for finite element simulations is here adopted in order to identify the stress fields and the constitutive parameters of heterogeneous materials. The change from an Airy's type approach to a Finite-Element approach to the CEGM allows investigating enhanced boundary conditions and more complex geometries.
In this application, we introduce the elastoplastic secant stiffness tensor B s . For a Prager's linear kinematic hardening model, the tensor B s is directly expressed as a function of the material properties (Young's modulus E, Poisson's ratio v for isotropic elasticity and shear modulus G for cubic elasticity, yield stress σ 0 and hardening coefficient h) and of the loading history. In the following, we briefly present the identification procedure, we illustrate its performance on various simulated data obtained under small perturbations and plane stress assumptions. Finally, we check its robustness with respect to measurement noise.

Identification procedure
Since we have interest in the sequel for thin and flat samples observed via in-plane DIC techniques, we focus on the identification of elastoplastic constitutive laws in a 2D framework (plane stress). In this context, a maximum of three parameters can be locally identified because we have only access to three local in-plane strain measurements related to three in-plane stress components.
The CEGM is based on the minimization of an energetic functional depending on two sets of parameters: the stress field and the mechanical material properties. This procedure can be applied to any identification problem and can be used both with data extracted from numerical simulations and with experimental measurements.
For a sequence of successive load steps (subscript 1 ≤ n ≤ N for each step), we denote by − → u m n the measured displacement field on a given region of interest and we consider an elastoplastic body governed by the set of Eqs. (1)(2)(3)(4): where σ c n represents a statically admissible stress field associated with the displacement − → u c n via the fourth order secant elastoplastic tensor B s n (corresponding to the Hooke tensor B e for an elastic step) for each load step n. It is worth noticing that for a heterogeneous material, all these quantities depend on the position.
The overall forces R j n are known for each time step n on the boundary ∂ j of . The free boundaries ∂ i satisfy the relations: On some boundaries ∂ u , we also impose that the average displacement − → u c n is equal to the average displacement − → u m n which not only allows to eliminate the rigid body motion but also to further constrain the identification problem and to reduce the influence of measurement noise.
The energetic functional can be expressed in its simplest form (small strain hypothesis, equilibrium): [1,N ] , B s n n∈ [1,N ] Here N is the overall number of time steps used for the identification and − → u c n a displacement field compatible with the equilibrium of the studied domain at load step n.
The identification procedure consists in minimizing E rc with respect to its two arguments − → u c n and B s n . The method can be applied to any behaviour for which an expression of the secant tensor B s n is available. Consequently, it can be used on reversible behaviours (linear and non-linear elasticity) or irreversible behaviours (viscoelasticity, elastoplasticity ...). When dealing with irreversible behaviours, it is necessary to take the loading history into account. In the case of elastoplasticity, this amounts to separate elastic loading steps from plastic ones. The method was numerically implemented to deal with elastoplastic behaviours and monotonic loadings. The consistency between the numerical implementation and the main hypotheses underlying the description of plastic flow (e.g. existence of a yield stress, isochoric plastic strain, normality rule) is ensured from the formulation, thus minimizing the set of parameters to identify.
In the case of cubic material, the Hooke tensor B e depends only on three elastic constants: E, v and G. According to [25], the fourth order secant elastoplastic tensor can be written at load step n as: where P is the constant mapping matrix: and γ n is the plastic multiplier that depends, for linear kinematic hardening, on the material parameters: with α + the positive part of a, α n the second invariant of the effective stress, X n the backstress tensor reached at the current load step: The expression of the secant modulus B s n plays a central role in the proposed method since it governs the definition of the plastic parameter K n which involves the plastic parameters to be identified. The secant tensor at time step n can also be expressed with respect to material parameters: Finally, since the plastic deformation at load step n is equal to ε p n = ε − B e−1 : σ c n , the plastic parameter K n , can be expressed as a function of two material parameters a and b: with a = 1 2h and b = σ 0 h in the case of a linear kinematic hardening. Note that when the load step n is purely elastic, the plastic strain is vanishing and the plastic parameter K n is equal to zero: the plastic secant tensor B s n is thus equal to the elastic tensor B e . Furthermore, the Eqs. (10) and (11) show that the plastic secant tensor depends on set of the elastoplastic parameters p = [E, v, G, σ 0 , h] and also on the norm of the plastic deformation ε p n so: where q are the phases (material domains) of the specimen and q p the material parameters of phase q.
To compute the ERC functional, several fields are to be defined: (1) the phase distribution (related to the material heterogeneity), (2) the stress fields (related to the development of structural effects in the specimen), and finally (3) the experimental displacement fields obtained by DIC.
These three fields are discretized using different meshes with adapted mesh size and shape functions. The meshes are "nested", the coarser being the material mesh, the finer being the DIC mesh, and the intermediate being the stress mesh. These three meshes are described with different shape functions. The stresses are determined through a FE computation using bilinear displacement elements and bilinear shape functions are used to perform the local DIC computation. The continuity of the measured displacement field is enforced by averaging the displacement vector on the mesh vertices and the mechanical properties are constant on each material domain. Moreover, it is possible to use different stress meshes in both elastic and plastic identification.
The 'plastic' mesh was a subdivision of the 'elastic' mesh in order to reduce the influence of the noise on the identification while maintaining a convenient description of the stress gradients. As conform meshes, these nested meshes simplify the transfers of fields from one mesh to another.
Finally, the direct simulations are performed by imposing constant vertical displacement on the upper boundary and blocking vertical displacement on the lower boundary. The left and right boundaries are stress free. Moreover we set the displacement of one point of the lower boundary to zero in order to suppress the possible horizontal rigid body motion.

Numerical method
Due to the convexity of the E rc functional, the minimization is performed in two steps, leading to the relaxation algorithm presented in Fig. 1a: the function is minimized with respect to the displacement field − → u c n associated with a statically admissible stress field, and then with respect to the secant tensor B s n . As already mentioned, we focus on a J2 elastoplastic model with kinematic hardening and no more than three material constants can be locally identified at each load step. The identification algorithm involves three steps (see Fig. 1b): (i) an elastic identification, (ii) a plasticity detection and (iii) a plastic identification.
The elastic and plastic parameters are thus identified separately and are based on the minimization of the E rc cost function. Nevertheless, in either situation, the first minimization is identical: it consists in a classical direct finite element computation for a known heterogeneous material under given boundary conditions. This minimization will thus not be discussed in the following. The elastic constants (i) are identified by minimizing the functional E rc with respect to the elastic tensor B e . The iterative minimization process starts with an initial set of parameters B 0 chosen by the user at n = 0 and taken equal to the elastic parameters identified on the previous load step when n > 1. The procedure is stopped using a convergence criterion defined on the norm of the correction in secant tensor. The plastic identification problem consists in determining the elastoplastic parameters involved in the secant stiffness tensor B s n . The plastic identification (step (iii)) is less direct since the secant tensor expression requires knowledge of the stress state. For all the plastic load steps, the secant tensor is initialized using the results of an elastic estimation B s 0,n , and gives a statically admissible stress state.
The plastic identification consists in minimizing E rc with respect to the plastic parameters (a and b for a linear kinematic hardening). Once the procedure has converged, we get the stress field σ c n and the optimal values (a opt , b opt ) that are directly related to the material parameters σ 0 and h. In this case, the minimization is performed numerically, due to the lack of analytical solution. For a given behaviour (elasticity or plasticity), and a given load step n, the identification algorithm presented in Fig. 1a is performed until convergence. In both cases, the convergence criterion is computed on the secant tensor The plasticity detection (step (ii)) consists in comparing the secant tensor B s N identified at the current load step N and the one identified at the previous (N − 1) load step B s N −1 : with b about 0.05. If this criterion is validated, we assume that the plasticity has started developing during the current load step N . The last elastic load step is denoted by N e . The plastic identification (step (iii)) is performed only for the load steps greater than N e (the overall loading is supposed to be monotonous).
Once the procedure has been completed, we get the statically admissible stress field for all the load steps and the set of identified material properties.

Validations
In this section, the efficiency of the proposed procedure is examined using reference simulated measurements obtained with the finite element code Comsol Multiphysics. Only simulated measurements are considered in this paper in order to focus only on the identification procedure performance and on the enhanced features of the formulation.
The in-plane components of the displacement field are extracted at the nodes of the finite element mesh and the global load levels are extracted on the outer boundaries. We have used different meshes for the direct computation and for the definition of the "measurement grid".
The errors on the identified parameters can have different origins. They can be related to errors on the measured displacement (i.e. DIC errors), on the geometry description (domain and boundaries), on the behavior law (description of the secant tensor B s n ), on the phase description (material mesh), or on the stress description (stress mesh and boundary conditions used for the stress computation). In this paper we address the influence of DIC errors, of phase description errors and, in a lesser extent of stress description errors. Since conform meshes are used, phase description errors are related both to stress description errors and to measured displacement errors.
The influence of the measured displacement field characteristics (mesh size and noise level) on the identification results is illustrated on numerical examples corresponding to homogeneous and heterogeneous materials subjected to a tensile test. It is well known that several error regimes can be encountered using DIC, noticeably the shape function mismatch regime and the ultimate error regime [26]. The shape function error is preponderant when the shape function used in the DIC formulation does not match the actual displacement field. It is governed by the first neglected term in the shape function and it tends to increase when the subset size is increased. The ultimate error is encountered when the shape function is rich enough to describe the real displacement field. It is governed by the image noise σ n , the average image gradient and the subset size d and it tends to decrease when the subset size increases [27]. These two aspects were separately investigated. The influence of the shape function mismatch error is examined using different displacement meshes leading to different mismatch error levels. The influence of the ultimate error is investigated solely by randomly perturbing the measured displacement field.

Elastic identification
The first test (specimen 1) is performed on an elastic bi-material sample: a soft circular isotropic inclusion (Young's modulus 100 GPa, Poisson ratio 0.15) is embedded in a stiff isotropic matrix (Young's modulus 210 GPa, Poisson ratio 0.3). Two types of identification mesh are used (Fig. 2): • An identification mesh perfectly consistent with the material domains (two identification domains D 1 and D 2 corresponding respectively to the inclusion and the matrix); • An identification mesh that does not match the material heterogeneity by splitting it into 400 domains D j with j = 1-400.
Only one load level is used in the identification procedure. The typical computing time on a Z820 workstation (Intel Xeon 2.40 GHz) computer is about 5 mins for a measurement mesh involving 1448 linear elements. For the first identification, Fig. 3a shows a good prediction of the parameter sets with a relative error of about 1% on the Young's modulus. In this case, the shape function mismatch error is small since the meshes used for the direct computation (considered as the "ground truth") and for the identification are identical.   The maximum error is about 4% of the maximum stress (295 MPa) and is located in the zones of maximum stress gradient.
For the second identification (involving 400 non-consistent material domains), the position of the inclusion is very well identified (see Fig. 3b).
This result shows the ability of the method to identify heterogeneous elastic properties without any a priori knowledge of the phase distribution. The error increases with the stress and strain gradients where the shape function mismatch error is the most significant. Furthermore, it can be observed that the error on the identified Young's modulus is concentrated above and under the inclusion where the deformation energy is minimal. This error is also important on the material domains intersecting the actual boundary between the inclusion and the matrix where the method tends to average the elastic constants of the two phases. The computational time depends on the number of unknowns which increases here from 4 in the previous case up to 800 in the case of 400 non-consistent material domains. It goes up to 8 mins on the same computer for the elastic identification.

Elastoplastic identification
The second test (specimen 2) concerns a standard tensile test performed at constant velocity on an isotropic elastoplastic material (Young's modulus 210 GPa, Poisson ratio 0.3, yield stress 300 MPa and hardening modulus 1 GPa) (Fig. 5a). The material parameters are identified using data associated with 5 load steps (2 steps in the elastic domain, and 3 in the plastic domain) (Fig. 5b). Although the material is homogeneous, the identification is made on 4 material domains (Fig. 5a) in order to demonstrate the ability of the procedure to identify the elastoplastic properties in several domains.
Identified parameters obtained from each zone are collected in Table 1. As can be noticed, the identified parameter values are very close to the reference values and are very similar from one identification domain to another.
The reference ("measured") stress fields presented in Fig. 6a are obtained by solving the direct problem whereas the stress fields presented in Fig. 6b are obtained by the inverse method. We notice a close similarity between the distributions and the orders of magnitude for this stress component. The procedure converges in few iterations. The identification of the parameters and of the given stress fields is very accurate.
In this case, different meshes were used for the direct computation and the identification: the former is triangular and the latter is quadrangular (see Fig. 6a, b). The mesh used to describe the "measured" displacement field (not represented here) is quadrangular but similar in size to the one used for the direct simulation thus limiting the bias introduced in the description of the "measured" displacement. The mesh used to compute the stress field is coarser than the one used for the direct simulation thus introducing a shape function mismatch. Here, the fact that the identification results are very close to the reference values show that the shape function mismatches is too small to generate significant discrepancies in the identification. The typical computing time on the same computer is about 6 mins for a measurement mesh involving 1448 linear elements. This computing time depends on the number of the load step used for the elasto-plastic identification.

Sensitivity to the initial set of parameters
To assess the sensitivity of the identification results to initial values, different starting values of the parameters are selected for the procedure. The identification is performed on specimen 2 and we check the number of iterations required for the procedure to converge. Table 2 shows that the identified parameters are very stable with respect to the chosen initial values. As mentioned earlier, the initial set of parameters is only used for the elastic identification of the first load step. No ad hoc initiation is required for the identification of plastic parameters.

Sensitivity to experimental noise
The robustness of the CEGM approach with respect to noise is evaluated using a set of simulated displacement fields disturbed by a white noise at different levels. For this purpose, we perform an identification on a homogeneous isotropic material submitted to a tensile test. The reference FE-displacement fields are corrupted by a white Gaussian noise with the amplitude γ. The noise level is chosen to γ = 0.01 pixel while the maximum displacement is 1.5 pixels. This value was chosen to be consistent with our classical test configurations: image noise σ n ≈ 0.54 grey levels (obtained using a HR16070MFLGEC camera, and a 16-bits acquisition mode), coarse speckle (speckle dots with 3-pixels radius, corresponding to an average image gradient ∇I 2 ≈ 90), and 20-pixels subset size (d). This value is consistent with the value σ n / d ∇I 2 ≈ 0.003 given in [27].
It can be seen that identification of all parameters is stable in presence of noise. As expected, the elastic constants are more corrupted by the noise level as, for a fixed noise level, the signal to noise ratio is smaller for the elastic identification than for the plastic one. Furthermore, the Poisson ratio is more sensitive to measurement noise.
The results presented in Table 3 are obtained using a stress identification mesh equal to the mesh used to obtain the direct problem and also without filtering the noise. But in our minimization algorithm we have several types of meshes that can have equal or different sizes. In order to reduce the problem of measurement noise in the identification while maintaining a good description of the stress gradients, these meshes are not identical but they are imbricated i.e. the mesh of the plastic identification is a subdivision of the elastic meshes.

Conclusions
In the present work, we use full-field measurements and the constitutive equation gap method to identify the spatial distribution of a set of material parameters associated with a J2 elastoplastic behaviour. The identification approach is based on the minimization of the CEG energy norm and allows the identification of a set of unknown parameters in any chosen zone without any prior knowledge of the distribution of the spatial heterogeneities.
We validate this approach on different materials in different situations (heterogeneous elastic materials, homogeneous plastic material involving heterogeneous stress fields). These different numerical tests give results which are in good agreement with the imposed value obtained using the direct problem. The identification accuracy strongly depends on the accuracy of the input data (measured load and displacement) and of the representativeness of the model (geometry, boundary conditions, material heterogeneity, behaviour law, stress distribution ...). Direct numerical simulations are performed in order to get the data required for the identification in a "perfectly controlled" situation. The results of these simulations are considered as a "ground truth" (thus neglecting the simulation error). The sensitivity of the identification results with respect to different parameters are investigated by perturbing the numerical solution. To restrict the scope of the study, we have focused here on three contributions: the measured displacement error, the shape function mismatch in the stress distribution and the description of the material heterogeneity. We have verified that typical error levels associated with the ultimate error regime of the displacement measurement led to relative errors smaller than 5% on the identification results (the elastic parameters being more sensitive to noise due to smaller signal to noise ratio). As expected, the shape function mismatch error is larger in stress concentration areas, and the identification error is more important in stress concentration zones. Special attention should be taken to choose a material mesh consistent with the phase distribution and to define a stress mesh fine enough to catch the stress concentrations. Finally, we have verified that the identification results are not affected by the choice of the initial guess on the elastic parameters required to start the procedure.
Future works will focus on the use of experimental data and the improvement of the method to deal with respect to multi-linear problems.