Extracting Constitutive Mechanical Parameters in Linear Elasticity Using the Virtual Fields Method Within the Ordinary State-Based Peridynamic Framework

This paper implements the virtual fields method within the ordinary state-based peridynamic framework to identify material properties. The key equations derived in this approach are based on the principle of virtual works written under the ordinary state-based peridynamic formalism for two-dimensional isotropic linear elasticity. In-house codes including a minimization process have also been developed to implement the method. A three-point bending test and two independent virtual fields arbitrarily chosen are used as a case study throughout the paper. The numerical validation of the virtual fields method has been performed on the case study by simulating the displacement field by finite element analysis. This field has been used to extract the elastic material properties and compared them with those used as input in the finite element model, showing the robustness of the approach. Noise analysis and the effect of the missing displacement fields on the specimen’s edges to simulate digital image correlation limitations have also been studied in the numerical part. This work focuses on pre-damage properties to demonstrate the feasibility of the method and provides a new tool for using full-field measurements within peridynamics with a reduced calculation time as there is no need to compute the displacement field. Future works will deal with damage properties which is the main strength of peridynamics.

Domain of an elastic body in the reference configuration

Introduction
Peridynamics (PD) reformulates continuum mechanics balance laws [47] through the nonlocal integration of the interactions between material points. Each material point interacts with its neighbors within a finite δ-radius region, the so-called horizon. Contrary to classical continuum mechanics, the PD governing equations are directly expressed in terms of force and displacement between material points, instead of spatial derivatives. This makes PD a suitable framework to address discontinuities, like cracks occurring in damage process [48]. Moreover, crack initiation and growth are addressed with only a single critical damage parameter that is triggered when a threshold is reached. This threshold can for example be related to the well-known free energy release rate [26,49] or to the J-integral value [52]. Therefore, the cracks propagate autonomously within the PD framework and no predefined path, nor other external criteria are needed. Different fracture modes can therefore be accounted for [33]. Extended finite elements, cohesive zone models, and PD predictions are compared against experimental data [1] for crack initiation and growth. The authors found that PD yields the more natural crack path, including branching and micro-branching. Many authors have attempted to clarify the horizon concept [7,8] since this is a major PD parameter. In addition, well-known and accepted constitutive theories and meth-ods derived within classical continuum mechanics have yet to be transposed into the PD framework [6,20].
Full-field measurements can be used to obtain every material parameter from a single, heterogeneous, mechanical test [4]. Full displacement fields can be obtained with surface measurement techniques like digital image correlation (DIC). Finite element model updating method [11,39,44], the constitutive equation gap method [25], the equilibrium gap method [10], and the virtual fields method [43] have been developed within classical continuum mechanics to obtain material parameters from such full-field measurements. These methods were also described in terms of flexibility of use, full-field displacement data accuracy, noise sensibility, computation time, etc. [2].
To the best of the authors' knowledge, only a handful of publications focused on fullfield measurements and inverse methods to identify constitutive parameters within PD. An approach which combines DIC and PD was developed [56] to determine full-field displacements for problems involving discontinuities. The principle is to use DIC where it is accurate, i.e., in continuous zones, and to add PD calculations in discontinuous regions. Thus, the full-field displacement solution is the combination of the DIC and PD calculations. An inverse method within state-based PD and using a Tikhonov regularization was developed [57] to extract the material properties and it provides promising prospects for damage characterization. However, it requires to compute the displacement field from the PD model at each step, which can be time consuming. Also, the robustness of the method against the initial guesses feeded in the minimization objective function and the selection or the effect of the horizon on the convergence were not studied.
The purpose of this work is to derive an equivalence to the virtual fields method (VFM) within the PD framework to identify material properties. A key point of the VFM, when compared with the other methods, is that no computation of the displacement field during the optimization process is needed. Moreover, this work focuses on predamage properties to demonstrate the feasibility of the peridynamics-virtual fields method (PD-VFM).
The paper is organized as follows. Section 2 recalls background information about the VFM within classical continuum mechanics and the PD framework. Section 3 explains the VFM using the theorem of virtual works formulated within the state-based PD framework. Section 4 presents a three-point bending simulation which is the case study used in the sequel. Section 5 describes the numerical implementation of the PD-VFM. Finally, Section 6 provides a numerical validation in which the displacement field of a three-point bending test has been simulated by finite element analysis and used to extract the material properties which have been then compared with those used as model input. The following notation is adopted in the paper: • Light-faced letters (e.g., a, α, A) denote constants or scalar quantities; • Bold-faced roman letters (e.g., a, A) denote first-order tensors; • Bold-faced lowercase greek letters (e.g., α) denote second-order tensors; • Non-italic bold-faced capital roman letters (e.g., A) denote fourth-order tensors; • Underlined light-faced lowercase letters (e.g., a) denote scalar states; • Underlined bold-faced capital letters (e.g., A) denote vector states; • "⊗" is the tensor product and ":" is the doubly contracted product; • Index notation and Einstein summation convention are used, unless specified otherwise.

Fundamental Equations
Consider an elastic body within a domain Ω bounded by ∂Ω and governed by the following equations: The Equilibrium Equation, Assuming Long-range Forces Are Negligible where σ is the Cauchy stress tensor, n the outward unit normal vector, and T the prescribed load acting on the surfaces S f .

The Kinematic Compatibility Equation, Within the Infinitesimal Strain Theory
where ε is the strain tensor, u the displacement vector, and u the prescribed displacement acting on the surface S u . S f and S u are such as S f ∪ S u = ∂Ω and S f ∩ S u = ∅.

The Constitutive Equation
where C is the stiffness tensor that depends on k constitutive parameters θ = (θ 1 , · · · , θ k ). For instance, an isotropic material depends on θ = (K, G) where K is the bulk modulus and G the shear modulus.

Set of Admissible Fields
where S and C are respectively the sets of statically admissible stress field and of kinematically admissible displacement field. It is worth noting here that the τ and v fields are independent of each other and not related through any constitutive model.

Principle of Virtual Works
The principle of virtual works (PVW) can be expressed from Eq. 1a as [43] − where u (assumed continuous) and ε = 1 2 ∇ · u + ∇ t · u are respectively the virtual displacement and strain fields. u ∈ C(0) means the virtual field is kinematically admissible such as u = 0 over S u in order to vanish the contribution of the unknown prescribed loading over S u . Both u and ε are independent of σ . Equation 6 can be divided into the following "physical" quantities: Thus, the PVW applied to a body in mechanical equilibrium yields It should be noted that the virtual quantity u , also known as virtual displacement, was implicitly assumed to have a length dimension. This provides the convenient "physical" interpretation of the principle of the virtual works for readers with engineering backgrounds. However, this virtual displacement has nothing to do with the actual displacement as it is a purely mathematical function bearing no relation whatsoever to physical quantities. Thus, the physical dimension of the different terms of the PVW depends on the physical dimensions of the virtual quantity.

VFM and Constitutive Theory Identification
The VFM is based on the PVW in which the strain field, denoted by ε, is experimentally determined using full-field measurement techniques. The nature of the material and its constitutive equation are assumed to be known. For example, for an elastic material where σ = C(θ 1 , · · · , θ k ) : ε, the governing equation becomes Then, at least k independent virtual fields must be chosen, among an infinite number of possibilities, to devise a system of k equations whose solution yields the k unknown constitutive parameters. The choice of these virtual fields is a key step. Three main approaches have been studied so far [2]: 1. Manually using a polynomial function (see the example below); 2. Automatically with special fields [30] that directly supply the constitutive parameters from the virtual work of the external loads. These special fields also reduce the technique's noise sensitivity [3]; 3. Piece-wise within the material [55]. Continuous lower-degree polynomial functions between each sub-region can be used. These latter avoid to use a high-degree polynomial defined on the whole geometry, and which would magnify the noise negative effects.
It should be noted that the volume integral in Eq. 8 requires the knowledge of the strain field in the solid. This is a drawback of the VFM since most of full-field measurement techniques, e.g., digital image correlation (DIC), provide in-plane strains over the specimen external surface. Plane stress, plane strain, or bending in thin plates loads are typically used to circumvent this limitation.
The VFM was applied to identify the material properties such as: -The in-plane linearly elastic properties of orthotropic composite plates from heterogeneous strain fields [29,31,55]; -The properties of a non-linearly elastic material assuming a polynomial stress/strain relationship [31]; -The bending rigidities of anisotropic thin plates within the framework of the Kirchhoff-Love theory [28,32]; -The stiffness and damping properties of thin viscoelastic isotropic vibrating plates [27]; -The properties of a material within a plastic behavior using an iterative procedure based on a minimization process as there is no direct relation between σ and the measured strains ε [55]; -The damage parameters [43] within the continuum damage mechanics [35].
All of these examples show that the VFM can be adapted to several kinds of problems. However, in all of them, the success of the VFM relies on the optimal choice of virtual fields.

Fundamental Equations
In peridynamics (PD) [6,47,52], a material point located at x interacts with all other material points within a region H x , centered on x and of a radius δ, which is referred to as the horizon (see Fig. 1). These material points in interaction with x are defined in the reference configuration Ω 0 and are not updated over time.
Let a material point be located at be the dual force density that a point at x exerts on a point at x. f v is a force per unit volume squared and has the following properties 's anti-symmetry results from the reciprocity principle. Thus, f v can be expressed using the bond force density t v as

Fig. 1 A material point located at
x interacts with all other material points within a region H = H x ∩ Ω 0 , centered on x and of a radius δ, which is referred to as the horizon [16] When assuming that long-range forces (e.g., gravity) are negligible, and similarly to Eqs. 1a and 2a, the PD equilibrium equation is F and u are the prescribed non-local boundary conditions acting on the boundary volume regions V f and V u . These are thin real material layers on the boundary of Ω and are such that is the external load applied as a body force density, i.e., b = 0 in Ω V f . This is different from classical continuum mechanics in which the boundary conditions are prescribed over the surface ∂Ω. Figure 1 shows that the locality of PD depends on the horizon radius δ. In elasticity, it is demonstrated that PD converges to classical continuum mechanics when the horizon approaches zero [51,52].

State-Based Theory
A mathematical object called states was developed [50,52] for the state-based theory. The states and their properties are presented in details in App. A of [20]. The main definitions are recalled here.
A vector state is a general mathematical object that maps vectors onto vectors. The wellknown second-order tensor can be viewed as a special vector state. In classical continuum mechanics, the constitutive model is expressed through a relationship between the Cauchy stress tensor σ and the strain tensor ε. Similarly, a PD constitutive model is expressed with a force state T linked to a deformation state Y.
Let x ∈ Ω 0 and x ∈ H x denote two material points in the reference configuration. Let ξ be the vector defined by ξ = x − x and called the bond connected to x. The family H of x is defined by where H and H x are schematized in Fig. 1. Let Y be the deformation state that transforms bonds ξ connected to x to their deformed images as where y(x) and y(x ) are, respectively, the deformed positions of x and x . Figure 2 illustrates these variables. It is further assumed that which means that two distinct points in the reference configuration Ω 0 are also distinct in the deformed configuration Ω t . Thus, the deformed direction vector state M is defined by The force state T is related to the bond force density t v through

Fig. 2
The deformation state Y transforms bonds ξ ∈ H to their deformed images. y(x) and y(x ) are, respectively, the deformed positions of x and x [14] Combining Eqs. 9a, 11a, and 16 leads to the PD equilibrium equation such as, ∀ x ∈ Ω 0 , App. B of [20] presents the three PD constitutive models which are the bond-based, the ordinary state-based and non-ordinary state-based theories. In the sequel, only the ordinary state-based approach is used since isotropic linearly elastic materials are considered.

2D Ordinary State-Based Model for Linearly Isotropic Elasticity
This section provides the constitutive models for an linearly isotropic elastic material under plane stress and plane strain (see the Appendix for details).

2D Plane Stress Constitutive Model
where α s and α d are the spherical and deviatoric PD parameters, e s and e d the spherical and deviatoric extension state, ν the Poisson's ratio, and q the PD weighted volume (see Eq. 37 and [50] for more details).

Discretization
Several approaches have been used to discretize the PD equilibrium Eq. 17a such as Galerkin finite element methods [9], Gauss quadrature [58], the spatial discretization [24,42], or the EMU (name of the first PD software implemented [40]) nodal discretization (ND). The last method (EMU-ND) consists of a midpoint (or one-point) quadrature in a Lagrangian spatial discretization. The EMU-ND has been chosen for its efficient load distribution scheme and acceptable computation time. Figure 3 shows the mesh of a reference position in 2D with regularly spaced nodes where h is the fixed mesh size. Each node i, located at x i in the dual mesh, is associated to a fixed volume V i . Assume that the volumes surrounding the nodes are non-overlapping (V i ∩ V j = ∅), and are recovering the reference volume V Ω , meaning that n i V i = V Ω . The discrete family H i of a node i is defined by: Thus, the discrete PD equilibrium equation under the EMU-ND scheme, and which has to numerically be solved for each node i, reads where b i is the external load applied to the node i. Recall that b i = 0 if node i is not on an external layer of the body where the external load is applied. Fig. 3 2D EMU-ND and its mesh with a constant mesh size h. V i is the surrounding volume of the node i located at x i and H i its discrete family. The integration over H x from Eq. 11a is now replaced by the discrete sum fron Eq. 21a over H i [15] One drawback reducing the accuracy of this discretization is the treatment of nodes j whose V j is partially inside the horizon of a node i and partially outside (see Fig. 3). Several algorithms were evaluated [46] to improve the accuracy of the discretization using partial areas taking into account the fraction that is outside the horizon. The most common is the PA-PDLAMMPS algorithm [42] which has been implemented in the code [19] presented in Section 5.

Principle of Virtual Works
Since Eq. 17a is valid for any x ∈ Ω 0 , each term can be "multiplied" (using the Euclidean inner product) by any continuous vector function u such as Integrating Eq. 22 over the whole domain Ω 0 yields the PVW stated as Assume that the virtual quantity u is kinematically admissible, i.e., u ∈ C(0), then, the contribution of the unknown force over the boundary volume constraint region V u vanishes. Equation 23 can be divided into the following "physical" quantities: where W int is the virtual work done by the internal dual force densities and W ext the virtual work done by the external forces. As for CCM version of the PVW, W int + W ext = 0.

Virtual Fields Method and Material Identification
The PD ordinary state-based constitutive equations for linearly isotropic elasticity depend non-linearly on the material properties (Eqs. 18a and 19a). Because of this non-linear relationship, the k sought material parameters are solution of a non-linear system. Rather than explicitly solving this non-linear system, the non-linear set of at least k independent equations were solved by minimizing a residual function Φ. Since linearly isotropic elastic materials are considered, Φ requires at least two independent virtual fields u (k) to prop-erly identify the material parameters (K, G). Thus, the residual function Φ, defined as a normalized squared residual between W int and W ext , was Note that the normalization was used only for comparison purpose. Also, it was assumed

The Virtual Displacement Fields Definition
Two trivial independent continuous virtual displacement fields u (1) and u (2) were defined: u (2) Fig. 4 Three-point bending setup in which a beam, of dimensions L × W × t and simply supported on two supports spaced of S, is submitted to an external force F . The coordinate system (O, e 1 , e 2 ) origin is located at the specimen's bottom face center [12] These virtual fields, depicted in Fig. 5, are kinematically admissible, i.e., they vanish on the supports:

Application of the Virtual Field Method
The internal virtual works W int were approximated from Eqs. 21a and 24a as The external virtual works W ext were calculated from Eq. 24b and are equal to: Finally, K and G were extracted using the residual function Φ from Eq. 25 and the minimization process described in Section 5. Figure 6 shows the algorithm state chart diagram used for the minimization process to extract the material properties from the actual displacement field. The algorithm has in input  State chart diagram of the algorithm to extract the material properties from the actual displacement field. The non-linear black box optimization software NOMAD [38], based on the MADS algorithm, is used for the minimization process the initial guess for the material properties (K; G) and the actual displacement field u. A uniformly distributed random value in a search interval I sufficiently large to the included material's expected properties was chosen as the initial guess. The gathered steps in the box (see Fig. 6) are repeated until the residual Φ is smaller or equal to a fixed tolerance . The main steps are described below:

Algorithm
1. The dual force density f v is computed from Eq. 21a using the trial material properties and u. For this step, the in-house state-based PD code PeriPyDIC [19] is used with the 2D constitutive modeling described in Section 2.2.3). Note that instead of solving a direct problem in which the displacement field is computed from the boundary conditions, here the dual force density f v is evaluated from a provided displacement field.

The internal virtual works W int are computed from Eq. 30 with the in-house code
PeriPyVFM [18], with the two virtual fields given in Eqs. 27 and 28. Note that PeriPyVFM is based on the discretized version of the equations given in Section 3 and can be extended to as many virtual fields as required. 3. The residual Φ defined in Eq. 25 is computed within PeriPyVFM using W int and the external virtual works W ext calculated from Eq. 31a.
As long as the residual Φ is greater than the tolerance , the "black box" optimization solver NOMAD [38], which implements the mesh adaptive direct search (MADS) algorithm, provides the new guess for the material properties. Thus, the three steps inside the box are repeated until black box minimization finds suitable values for the material properties and the algorithm terminates. Table 1 summarizes the versions of the software used in this paper. PeriPyDIC 0.1 [19] PeriPyVFM 0.1 [18]

Choice of the Horizon δ
Choosing the ratio δ/h of horizon and nodal spacing to gain convergence is still an open question within the PD community [7,8]. Local limits and asymptotically compatible discretizations were considered in [23] and [22]. It is shown in [54] that, if the nodal spacing h → 0 decays faster than the horizon δ → 0, the EMU-ND converges to the continuum local partial differential equation (PDE) limit. An empirical choice of the horizon for a specific material is to adjust the horizon to fit experimental results. Convergence analysis can also be performed as described in [8] in which three types of convergence are presented: the δ-, m-, and (δm)-convergence where m = δ h . The approach used in this paper assumes the use of experimental data. In that case, the nodal spacing h is more or less prescribed by the measuring equipment resolution, e.g., the DIC system. As the horizon is such as δ = mh, m is the only parameter to adjust with respect to the residual function Φ. Thus, a three-point bending test simulated by finite element with a mesh size of h = 0.25 mm was used to calibrate m. The residual Φ was computed for the exact material properties for different value of m. Figure 7 plots Φ as a function of m. It can be seen that Φ is large for m ≥ 2, gets smaller for the next values, and is the smallest for m = 7. Then, Φ increases for the next two values. This kind of non-stabilized behavior has been shown in [21]. Therefore, m = 7 was used in the sequel to extract the material properties (K, G) using the algorithm presented in Fig. 6. Note that this is not a convergence study but only a simplistic way to choose the horizon δ.

Actual Displacements Simulated by Finite Element Analysis
The virtual three-point bending test presented in Fig. 4 was modelled in finite element analysis (FEA) with ANSYS . The setup size parameters were L = 128mm, W = 32 mm, t = 12.7 mm, and S = 96mm. The geometry was discretized using the PLANE182 2-D Structural Solid element with a mapped meshing of h = 0.25 mm. The element was used as a plane stress element and was defined by four nodes having two translations at each node in e 1 and e 2 directions. A null displacement condition was applied to the nodes located on the two supports and the beam was loaded with a force of magnitude F of 1 000 N applied to node A shown in Fig. 4. The Young's modulus E was of 4 000 MPa and the Poisson's ratio ν of 0.3, which yields a bulk modulus K = 3 333.33 MPa and a shear modulus G = 1 538.46 MPa. Figure 8 a and b respectively show the components u 1 and u 2 for the simulated displacement field u FEA .

Extracted Material Properties Using the Displacement Field u FEA
The displacement field u FEA was input in the algorithm described in Section 5 to identify the material properties. The search interval for the sought (K, G) was I = [1 000; 10 000] 2 . This interval includes the entered values K = 3 333.33 MPa and G = 1 538.46 MPa. Based on Eq. 31a, the external virtual works were W (1) ext = 48 000 N mm and W (2) ext = 0 N mm. The extracted material properties are presented in Table 2.
When the initial guesses for K and G were exactly in the mid of the search interval I , i.e., (K; G) = (5 000; 5 000), the relative errors with respect to the material properties fed in the FEM were K = 0.63% for the bulk modulus and G = −0.44% for the shear modulus. The extracted properties using the state-based PD-VFM are in good agreement with those utilized in the FEA.
The minimization was run 14 more times with different initial guesses distributed within the search interval to investigate the method's robustness and stability. Table 2 gathers the results and shows that all the extracted material properties converge to the same values (they only differ for the fourth decimal point), even though the number of evaluations was different. Therefore, the minimization of the residual Φ seems to be stable and means that exactly one global minimum must exist within the search interval.

Noise Analysis
The actual displacement field u FEA used to identify the properties was optimal with respect to bias and noise as this is a simulated field. Thus, the slight differences between the actual   The first column is the noise's standard deviation used. The second and third columns represent the extracted material properties after minimization. The fourth and fifth columns show the relative error with respect to the material properties use in the FEA, which are (K; G) = (3 333.33; 1 538. 46) material properties fed in the FEM and those extracted could be attributed to several sources such as the error induced with the discretization scheme, the so-called PD surface effect [36], the numerical round off errors or the virtual field chosen which can be more or less sensitive to the previous errors effects [43].
Real displacement fields are disturbed both by bias and noise [53]. Contrary to the DIC bias that can be reduced or eliminated with a proper setup and parameters (good calibration, no contamination or dust, no aliasing), the noise is unavoidable even though it can be minimized with a careful setup (good focus and speckle pattern, contrast, lighting, stereo-angle, lens selection). A noise analysis simulating experimental data was performed by adding a Gaussian noise term to disturb u FEA as where g NOISE is the disturbing displacement vector whose components g 1 and g 2 follow a Gaussian distribution of mean μ and standard deviation . It was assumed that μ = 0.0 mm, and lied between 0.0 and 25 × 10 −4 mm which is a typical standard deviation range measured with the DIC system used in the Laboratory for Multiscale Mechanics (LM2). Table 3 provides the extracted material properties after adding Gaussian noise. The relative error remained under 1.4% in absolute value. Therefore, the virtual fields chosen seem stable with respect to noise. This could be different with other virtual fields. Also, it should Fig. 9 Missing nodes all around the specimen's edges [13] Fig. 10 Influence of missing data: the black line shows the linear regression between the residual Φ with respect to the missing distance around the specimen be noted that the bulk modulus K is more sensitive to Gaussian noise than the shear modulus G. The conclusion from this noise analysis is similar that presented in [34], in which the VFM is applied within CCM.

Missing Data
Obtaining displacement fields on the exact specimen edges is hardly possible with DIC and the missing data due to this experimental limitation could have a significant influence on the identification results, depending the the virtual fields chosen [43]. Lines of elements that were placed of a multiple of h = 0.25 mm have been removed from the specimen's edges (see Fig. 9) before extracting the material properties to simulate the effect of the missing data. Figure 10 presents the effects of these missing nodes on the residual Φ. A linear relation can be highlighted, showing the necessity to obtain data as close as possible to the free edges The first column is the missing distance on the edge (see Fig. 9). The second and third columns represent the extracted material properties after minimization. The fourth and fifth columns show the relative error with respect to the material properties used in the FEA, which are (K; G) = (3 333.33; 1 538. 46) to accurately extract the material properties. Table 4 provides the error in the identification of K and G, showing that the relative errors quickly increase for the bulk modulus K. This means the identification of K is more sensitive to missing data than the shear modulus G using the three-point bending experimental setup and the virtual fields defined in Eqs. 27 and 28.

Conclusion
The main contributions of this study are as follows: 1. The VFM within PD has been derived to extract the elastic properties of a linearly isotropic material. Open-source python codes [18,19] including a minimization process have also developed to implement the PD-VFM. Contrary to most other approaches existing in the PD literature, one of the benefits of the PD-VFM is that no computation of the displacement field during the minimization process is required. 2. The PD-VFM has been tested on a three-point bending test using undisturbed displacement fields generated by FEA and two independent virtual fields arbitrarily chosen. The material properties extracted with the PD-VFM were in a quasi-perfect agreement with those utilized in the FEA. This numerical validation of the PD-VFM shows the robustness of the approach. Both the stability and the sensibility against noise and missing data have also been studied. The latter highlighted a linear relation between these missing data and the residual Φ, showing the necessity to obtain data as close as possible to the free edge to accurately extract the material properties.
Even though the numerical validation of the PD-VFM led to remarkable results, meaning the material properties identification within less than 0.75% error, there is still a scope for enhancing the solver performance and reduce the calculation time. Also, the authors tried to apply the PD-VFM to experimental DIC measurements performed on a polycarbonate beam under a three-point bending test. However, the residual Φ computed using the actual material properties (measured per ASTM D638-14) was too high and then did not allow to extract the material properties. This can be explained by the difficulty of reducing the missing data due to the DIC measurements and the high sensibility of the bulk modulus K with respect to this missing information and virtual fields chosen. Based on the full-field measurement literature [5,43], the authors would recommend to design a well-suited test configuration for the PD-VFM that would meets the following criteria: -A large heterogeneity of the displacement fields. This is the spirit of the full-field measurement methods, by opposition to standard mechanical tests that are based on homogeneous and averaged fields; -A good sensitivity of the displacement fields with respect to the material properties to identify.

A.1 Elastic Material
An elastic material is defined by its free energy density ψ which only depends on Y. Define ψ = Ψ (Y) as the strain energy density. Then, the force state T is defined by

A.2 Isotropy
A material is isotropic if and only if, for any proper orthogonal second-order tensor Q T(YQ) ξ = T(Y) Qξ .
Physically, Eq. 34 reflects the fact that the force state is invariant by applying rotations before deformations.

A.3 Isotropic Elastic Solid
Let ϑ be a non-local volume dilatation [36,37,45,52] such as w is the scalar influence state used to weight the bonds influence on the force state calculation and x = |X ξ | = ξ 2 . q is the weighted volume such as 2ν − 1 ν − 1 in 2D plane stress 2 in 2D plane strain (37) ν is the Poisson's ratio. Finally, e is the scalar extension state defined by The scalar extension state can be divided into a spherical part e s and a deviatoric part e d such as e = e s + e d = x 3 ϑ(e) + e d .
The spherical part represents the isotropic expansion while the deviatoric part depicts shear deformations.
Suppose an isotropic elastic material defined by its strain energy density Ψ Ψ (e) = 1 2 α s (w e s ) • e s + 1 2 α d (w e d ) • e d (40) where α s and α d are the PD parameters [20]. They are calibrated to the classical isotropic elastic material properties (K and G, the bulk and shear modulus) by equating the PD strain energy density Ψ to the classical strain energy density for an arbitrary homogeneous deformation under the infinitesimal strain theory [50]. One finds [20,36,37,45] that the PD parameters are: 1 q (8G) in 2D plane stress or plane strain (42) The calibration above also assumes that the material points have a complete neighborhood of size 2δ (i.e., H = H x ) filled with the same material. This assumption is not verified within a distance of δ next to a free surface or an interface with another materials and implies that Eqs. 41  represents an increment of the function of state resulting from an incremental change to the scalar extension state e. This incremental change to the scalar extension state is related to the incremental change to the deformation state Y as follows Thus, combining Eqs. 38, 33, 43a, and 44 leads to the constitutive model of an isotropic elastic material given by the force state Equation 45 shows that the constitutive model of an isotropic elastic material corresponds to an ordinary material since T = t