Extracting elastic constitutive parameters using the VFM within peridynamics

,


Introduction
Peridynamics (PD) reformulates continuum mechanics balance laws [Silling, ] through the non-local 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 [Silling, ].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 [Silling and Askari, , Foster et al., ] or to the J-integral value [Silling and Lehoucq, ].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 [Hu et al., ].Extended finite elements, cohesive zone models and PD predictions are compared against experimental data [Agwai et al., ] 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 [Bobaru et al., , Bobaru and Hu, ] since this is a major PD parameter.In addition, well-known and accepted constitutive theories and methods derived within classical continuum mechanics have yet to be transposed into the PD framework [Bobaru et al., , Delorme et al., ].
Full-field measurements can be used to obtain every material parameter from a single, heterogeneous, mechanical test [Avril and ] 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, fullfield displacement data accuracy, noise sensibility, computation time, etc. [Avril et al., ].
To the best of the authors' knowledge, only a handful of publications focused on full-field measurements and inverse methods to identify constitutive parameters within PD.An approach which combines DIC and PD was developed [Turner, ] 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 [Turner et al., ] 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 to the other methods, is that no computation of the displacement field during the optimization process is needed.Moreover, this work focuses on pre-damage properties to demonstrate the feasibility of the Peridynamics-Virtual Fields Method (PD-VFM).
The paper is organized as follows.Sec.recalls background information about the VFM within classical continuum mechanics and the PD framework.Sec.explains the VFM using the theorem of virtual works formulated within the state-based PD framework.Sec.presents a three-point bending simulation which is the case study used in the sequel.Sec.describes the numerical implementation of the PD-VFM.Finally, Sec.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 to 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.

Background . Virtual Fields Method (VFM) within classical continuum mechanics . . 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 (PVW)
The Principle of Virtual Works (PVW) can be expressed from Eq. ( ) as [Pierron and Grédiac, ] − 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 σ.Eq. ( ) 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 depend 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 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 [Avril et al., ]: . manually using a polynomial function (see the example below); . automatically with special fields [Grédiac et al., a] that directly supply the constitutive parameters from the virtual work of the external loads.These special fields also reduce the technique's noise sensitivity [Avril et al., ]; . piece-wise within the material [Toussaint et al., ].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. ( ) 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 [Grédiac et al., , Grédiac et al., b, Toussaint et al., ]; • the properties of a non-linearly elastic material assuming a polynomial stress/strain relationship [Grédiac et al., b]; • the bending rigidities of anisotropic thin plates within the framework of the Kirchhoff-Love theory [Grédiac, , Grédiac et al., ]; • the stiffness and damping properties of thin viscoelastic isotropic vibrating plates [Giraudeau and Pierron, ]; • 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 ε [Toussaint et al., ]; • the damage parameters [Pierron and Grédiac, ] within the continuum damage mechanics [Ladeveze and LeDantec, ].
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. .

Peridynamic (PD) framework . . Fundamental equations
In peridynamics [Silling, , Silling and Lehoucq, , Bobaru et al., ], 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. ).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 x in H x .Let f v = f v (x , x) 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 Eq. ( b)'s anti-symmetry results from the reciprocity principle.Thus, f v can be expressed using the bond force density t v as When assuming that long-range forces (e.g., gravity) are negligible, and similarly to Eqs. ( ) and ( ), the PD 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 [Delorme, e] 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 b 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 ∂Ω. . .

State-based theory
A mathematical object called states was developed [Silling et al., , Silling and Lehoucq, ] for the state-based theory.The states and their properties are presented in details in App.A of [Delorme et al., ].The main definitions are recalled here.
A vector state is a general mathematical object that maps vectors onto vectors.The well-known secondorder 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. .
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 .Fig. 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 Combining Eqs. ( ), ( a) and ( ) leads to the PD equilibrium equation such as, App. B of [Delorme et al., ] 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.

D 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 App.A for details).

D 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. ( ) and [Silling et al., ] for more details).

. . Discretization
Several approaches have been used to discretize the PD equilibrium equation ( ) such as Galerkin finite element methods [Chen and Gunzburger, ], Gauss quadrature [Weckner and Emmrich, ], the spatial discretization [Emmrich and Weckner, , Parks et al., ] or the EMU (name of the first PD software implemented [Littlewood, ]) 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.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 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.
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. ).Several algorithms were evaluated [Seleson, ] 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 [Parks et al., ] which has been implemented in the code [Delorme et al., ] presented in Sec. .

The Virtual Fields Method within Peridynamics . Principle of Virtual Works
Since Eq. ( ) 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. ( ) over the whole domain Ω 0 yields the PVW stated as Figure : D 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. ( ) is now replaced by the discrete sum fron Eq. ( ) over H i [Delorme, d].
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.Eq. ( ) 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.( ) and ( )).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 properly 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 that u (k) are

Case study: three-point bending test . Geometry and boundary conditions
Fig. describes a virtual three-point bending setup in which a beam, of dimensions L × W × t (length × width × thickness) and simply supported on two supports spaced of S. The coordinate system (O, e 1 , e 2 ) origin lies in the specimen's middle, on its bottom line.The beam is submitted to an external force applied to the specimen's top surface center as .

The virtual displacement fields definition
Two trivial independent continuous virtual displacement fields u (1) and u (2) were defined: These virtual fields, depicted in Fig. , 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. ( ) and ( a) as The external virtual works W ext were calculated from Eq. ( b) and are equal to: Finally, K and G were extracted using the residual function Φ from Eq. ( ) and the minimization process described in Sec. .origin is located at the specimen's bottom face center [Delorme, f]. .The dual force density f v is computed from Eq. ( ) using the trial material properties and u.For this step, the in-house state-based PD code PeriPyDIC [Delorme et al., ] is used with the D constitutive modeling described in Sec. . .).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.

Minimization process implementation . Algorithm
. The internal virtual works W int are computed from Eq. ( ) with the in-house code PeriPyVFM [Delorme and Diehl, ], with the two virtual fields given in Eqs. ( ) and ( ).Note that PeriPyVFM is based on the discretized version of the equations given in Sec. and can be extended to as many virtual fields as required.
. The residual Φ defined in Eq. ( ) is computed within PeriPyVFM using W int and the external virtual works W ext calculated from Eq. ( ).
As long as the residual Φ is greater than the tolerance , the "black box" optimization solver NOMAD [Le Digabel, ], 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.  .

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 [Bobaru et al., , Bobaru and Hu, ].Local limits and asymptotically compatible discretizations were considered in [Du and Tian, ] and [Du, ].It is shown in [Tian and Du, ] 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 [Bobaru et al., ] 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 threepoint 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.Fig. 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 [Diehl et al., ].Therefore, m = 7 was used in the sequel to extract the material properties (K, G) using the algorithm presented in Fig. .Note that this is not a convergence study but only a simplistic way to choose the horizon δ.

Numerical validation . Actual displacements simulated by Finite Element Analysis
The virtual three-point bending test presented in Fig. was modelled in finite element analysis (FEA) with ANSYS ® .The setup size parameters were: L = 128 mm, W = 32 mm, t = 12.7 mm and S = 96 mm.The geometry was discretized using the PLANE -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 .

Extracted material properties using the displacement field
The displacement field u FEA was input in the algorithm described in Sec. 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.33MPa and G = 1 538.46MPa.Based on Eq. ( ), the external virtual works were W   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 statebased PD-VFM are in good agreement with those utilized in the FEA.
The minimization was run more times with different initial guesses distributed within the search interval to investigate the method's robustness and stability.Table 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 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 [Le and Bobaru, ], the numerical round off errors or the virtual field chosen which can be more or less sensitive to the previous errors effects [ Pierron and Grédiac, ].
Real displacement fields are disturbed both by bias and noise [Sutton et al., ].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 2.5 × 10 −4 mm which is a typical standard deviation range measured with the DIC-system used in the Laboratory for Multiscale Mechanics (LM ).
Table : Extracted material properties from u FEA .15 initial guesses (K; G) ∈ [1 000; 10 000] 2 were used and provided in the second and third columns.The fourth and fifth columns represent the extracted material properties after minimization.The sixth and seventh 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).The last column shows the number of evaluations performed with the NOMAD solver.

Initial guesses Minimization results
Relative error .

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 [Pierron and Grédiac, ].Lines of elements that were placed of a multiple of h = 0.25 mm have been removed from the specimen's edges (see Fig. ) before extracting the material properties to simulate the effect of the missing data.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. ( ) and ( ).

Conclusion
The main contributions of this study are as follows: . The VFM within PD has been derived to extract the elastic properties of a linearly isotropic material.
Open-source python codes [Delorme and Diehl, , Delorme et al., ] 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.
. The PD-VFM has been tested on a three-point bending test using undisturbed displacement fields Table : Error in the identification of K and G due to Gaussian noise.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

Figure :
Missing nodes all around the specimen's edges [Delorme, b].
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 .% 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 D -) 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 [Belhabib et al., ,Pierron and Grédiac, ], the authors would recommend to design a well-suited test configuration for the PD-VFM that would meets the following criteria: Table : Error in the identification of K and G due to the missing data on the specimen free edges.The first column is the missing distance on the edge (see Fig. ).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

Ωσ
: ε dV Virtual work done by the internal forces W ext = S T • u dS Virtual work done by the external forces.

Figure
Figure :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[Delorme,  e]

Fig. shows
Fig.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[Silling and Lehoucq,  , Silling and Lehoucq,  ].

Figure :
Figure :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[Delorme,  c].

Figure :
Figure :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[Delorme,  a].
Shape of the second virtual field u(2)

Figure :
Figure :Virtual displacements applied to the three-point bending test.The coordinate system (O, e 1 , e 2 ) origin is located at the specimen's bottom face center[Delorme,  f].

Fig. shows the
Fig. 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 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 include material's expected properties was chosen as the initial guess.The gathered steps in the box (see Fig. ) are repeated until the residual Φ is smaller or equal to a fixed tolerance .The main steps are described below:

Start
Figure :State chart diagram of the algorithm to extract the material properties from the actual displacement field.The non-linear black box optimization software NOMAD[Le Digabel,  ], based on the MADS algorithm, is used for the minimization process.

Figure :
Figure :Residual Φ for the different m-values with h = 0.25 mm.The objective of this simplistic analysis is to find the best m-value to use in the PD model in order to obtain the smallest residual Φ for the prescribed mesh-size.
Fig. .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.33MPa and a shear modulus G = 1 538.46MPa.Figs. a and b respectively show the components u 1 and u 2 for the simulated displacement field u FEA .
000 N mm and W (2) ext = 0 N mm.The extracted material properties are presented in Table .
u 1 -component of the displacement field u FEA .(b) u 2 -component of the displacement field u FEA .

Figure :
Figure : Components u 1 (a) and u 2 (b) of the simulated displacement field u FEA .

Fig
Fig. 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 to accurately extract the material properties.Table 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. ( ) and ( ).

Figure :
Figure :Influence of missing data: the black line shows the linear regression between the residual Φ with respect to the missing distance around the specimen.
Pierron,  ].Full displacement fields can be obtained with surface measurements techniques like Digital Image Correlation (DIC).Finite Element Model Updating Method[Cottin et al., , Rouger et al., , Lecompte et al., ], the Constitutive Equation Gap Method [Florentin and Lubineau, ], the Equilibrium Gap Method [Claire et al., ] or the Virtual Fields Method [Pierron and Grédiac, Table summarizes the versions of the software used in this paper.

Table :
Software used for the implementation of the state-based PD-VFM.
Table 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 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[Kramer and Scherzinger,  ], in which the VFM is applied within CCM.