This paper implements the Virtual Fields Method within the ordinary state based peridynamic frame- work to identify material properties. The key equations derived in this approach are based on the princi- ple 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 /uniFB01elds arbitrarily chosen are used as a case study throughout the paper. The numerical validation of the virtual /uniFB01elds method has been performed on the case study by simulating the displacement /uniFB01eld by /uniFB01nite element analysis. This /uniFB01eld has been used to extract the elastic material properties and compared them to those used as input in the /uniFB01nite element model, showing the robustness of the approach. Noise analysis and the e/uniFB00ect of the missing displacement /uniFB01elds 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-/uniFB01eld measurements within peridynamics with a reduced calculation time as there is no need to compute the displacement /uniFB01eld. Future works will deal with damage properties which is the main strength of peridynamics. 1 Nomenclature Material properties E Young’s modulus [ MPa ]  Poisson’s ratio K Bulk modulus [ MPa ] G Shear modulus [ MPa ] Geometrical parameters h Mesh size / Nodal spacing [ mm] L Length of the beam [ mm] S Distance between the pins to supporting pins [ mm] Sf Surface with prescribed loading [ mm2] Su Surface with prescribed displacement [ mm2] t Thickness [ mm] Vf Volume with prescribed loading [ mm3] Vu Volume with prescribed displacement [ mm3] W Width of the beam [ mm] Classical continuum mechanics tensors  Cauchy stress tensor [Components in MPa ] " In/uniFB01nitesimal strain tensor b" Actual in/uniFB01nitesimal strain tensor "?Virtual in/uniFB01nitesimal strain tensor C Sti/uniFB00ness tensor [Components in MPa ] Peridynamic objects  Peridynamic horizon [ mm] H Family of a material point Hi Discrete family of a node i Hx Spherical region centered on x q Peridynamic weighted volume  Peridynamic bond fv Dual force density [Components in N mm6] tv Bond force density [Components in N mm6] M Deformed direction state T Force state [Components in N mm6] Y Deformation state [Components in mm] w In/uniFB02uence state esSpherical extension state [Magnitude in mm] edDeviatoric extension state [Magnitude in mm] sSpherical peridynamic parameter dDeviatoric peridynamic parameter /two.pnum Position and displacement vectors x Reference position vector [Components in mm] y Deformed position vector [Components in mm] u Displacement vector [Components in mm] u Prescribed displacement vector [Components in mm] bu Actual displacement vector [Components in mm] u?Virtual displacement vector [Components in mm] Force vectors b Body force density vector [Components in N mm3] F Force vector [Components in N] F Prescribed force vector [Components in N] T Prescribed surface force vector [Components in N mm2] Other Symbols I Search interval C Set of kinematically admissible displacement /uniFB01elds S Set of statically admissible displacement /uniFB01elds W? int Internal virtual work within classical continuum mechanics [ N mm ] W? ext External virtual work within classical continuum mechanics [ N mm ] W? int Internal virtual work within peridynamics [ N mm ] W? ext External virtual work within peridynamics [ N mm ]  Tolerance to end the minimization process K Relative error related on the bulk modulus G Relative error related on the shear modulus  Mean of a Gaussian distribution  Residual function to minimize in the virtual /uniFB01eld method  Standard deviation of a Gaussian distribution Domain of an elastic body 0 Domain of an elastic body in the reference con/uniFB01guration /three.pnum /one.pnum Introduction Peridynamics (PD) reformulates continuum mechanics balance laws [Silling, /two.pnum/zero.pnum/zero.pnum/zero.pnum] through the non-local integration of the interactions between material points. Each material point interacts with its neighbors within a /uniFB01nite -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, /two.pnum/zero.pnum/one.pnum/six.pnum]. 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, /two.pnum/zero.pnum/zero.pnum/five.pnum,Foster et al., /two.pnum/zero.pnum/one.pnum/one.pnum] or to the J-integral value [Silling and Lehoucq, /two.pnum/zero.pnum/one.pnum/zero.pnum]. Therefore, the cracks propagate autonomously within the PD framework and no prede/uniFB01ned path, nor other external criteria are needed. Di/uniFB00erent fracture modes can therefore be accounted for [Hu et al., /two.pnum/zero.pnum/one.pnum/five.pnum]. Extended /uniFB01nite elements, cohesive zone models and PD predictions are compared against experimental data [Agwai et al., /two.pnum/zero.pnum/one.pnum/zero.pnum] 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., /two.pnum/zero.pnum/zero.pnum/nine.pnum, Bobaru and Hu, /two.pnum/zero.pnum/one.pnum/two.pnum] 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., /two.pnum/zero.pnum/one.pnum/six.pnum, Delorme et al., /two.pnum/zero.pnum/one.pnum/seven.pnum]. Full-/uniFB01eld measurements can be used to obtain every material parameter from a single, heterogeneous, mechanical test [Avril and Pierron, /two.pnum/zero.pnum/zero.pnum/seven.pnum]. Full displacement /uniFB01elds can be obtained with surface measure- ments techniques like Digital Image Correlation (DIC). Finite Element Model Updating Method [Cottin et al., /one.pnum/nine.pnum/eight.pnum/four.pnum, Rouger et al., /one.pnum/nine.pnum/nine.pnum/one.pnum, Lecompte et al., /two.pnum/zero.pnum/zero.pnum/seven.pnum], the Constitutive Equation Gap Method [Florentin and Lu- bineau, /two.pnum/zero.pnum/one.pnum/zero.pnum], the Equilibrium Gap Method [Claire et al., /two.pnum/zero.pnum/zero.pnum/four.pnum] or the Virtual Fields Method [Pierron and Grédiac, /two.pnum/zero.pnum/one.pnum/two.pnum] have been developed within classical continuum mechanics to obtain material parameters from such full-/uniFB01eld measurements. These methods were also described in terms of /uniFB02exibility of use, full- /uniFB01eld displacement data accuracy, noise sensibility, computation time, etc. [Avril et al., /two.pnum/zero.pnum/zero.pnum/eight.pnum]. To the best of the authors’ knowledge, only a handful of publications focused on full-/uniFB01eld measurements and inverse methods to identify constitutive parameters within PD. An approach which combines DIC and PD was developed [Turner, /two.pnum/zero.pnum/one.pnum/four.pnum] to determine full-/uniFB01eld displacements for problems involving discontinu- ities. 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-/uniFB01eld displacement solution is the combination of the DIC and PD calculations. An inverse method within state-based PD and using a Tikhonov regularization was devel- oped [Turner et al., /two.pnum/zero.pnum/one.pnum/five.pnum] to extract the material properties and it provides promising prospects for damage characterization. However, it requires to compute the displacement /uniFB01eld 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 e/uniFB00ect 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 /uniFB01eld 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. /two.pnum recalls background information about the VFM within classical continuum mechanics and the PD framework. Sec. /three.pnum explains the VFM using the theorem of virtual works formulated within the state-based PD framework. Sec. /four.pnum presents a three-point bending simulation which is the case study used in the sequel. Sec. /five.pnum describes the numerical implementation of the PD-VFM. Finally, Sec. /six.pnum provides a numerical validation in which the displacement /uniFB01eld of a three-point bending test has been simulated by /uniFB01nite 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 /uniFB01rst-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; /four.pnum 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 speci/uniFB01ed otherwise. /two.pnum Background /two.pnum./one.pnum Virtual Fields Method (VFM) within classical continuum mechanics /two.pnum./one.pnum./one.pnum 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 (r=0in (/one.pnuma) n=TonSf@ (/one.pnumb) whereis the Cauchy stress tensor, nthe outward unit normal vector and Tthe prescribed load acting on the surfaces Sf. The kinematic compatibility equation, within the in/uniFB01nitesimal strain theory 8 >>< >>:"=1 2 ru+rtu in (/two.pnuma) u=uonSu@ (/two.pnumb) where"is the strain tensor, uthe displacement vector and uthe prescribed displacement acting on the surfaceSu.SfandSuare such asSf[Su=@ andSf\Su=;. The constitutive equation (=C:"in (/three.pnuma) C=C() (/three.pnumb) where Cis the sti/uniFB00ness tensor that depends on kconstitutive parameters = (1;;k). For in- stance, an isotropic material depends on = (K;G )whereKis the bulk modulus and Gthe shear modulus. Set of admissible /uniFB01elds S(T) = jr =0in and n=TonSf (/four.pnum) C(u) =fvjv=uonSug. (/five.pnum) whereSandCare respectively the sets of statically admissible stress /uniFB01eld and of kinematically admissible displacement /uniFB01eld . It is worth noting here that the andv/uniFB01elds are independent of each other and not related through any constitutive model. /five.pnum /two.pnum./one.pnum./two.pnum Principle of Virtual Works (PVW) The Principle of Virtual Works (PVW) can be expressed from Eq. (/one.pnum) as [Pierron and Grédiac, /two.pnum/zero.pnum/one.pnum/two.pnum] Z :"?dV+Z STu?dS= 08u?2C(0) (/six.pnum) whereu?(assumed continuous) and "?=1 2(ru?+rtu?)are respectively the virtual displacement and strain /uniFB01elds. u?2 C(0)means the virtual /uniFB01eld is kinematically admissible such as u?=0over Suin order to vanish the contribution of the unknown prescribed loading over Su. Bothu?and"?are independent of . Eq. (/six.pnum) can be divided into the following “physical” quantities: 8 >>>>< >>>>:W? int=Z :"?dV Virtual work done by the internal forces W? ext=Z STu?dS Virtual work done by the external forces. Thus, the PVW applied to a body in mechanical equilibrium yields W? int+W? ext= 08u?2C(0). (/seven.pnum) 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 di/uniFB00erent terms of the PVW depend on the physical dimensions of the virtual quantity. /two.pnum./one.pnum./three.pnum VFM and constitutive theory identi/uniFB01cation The VFM is based on the PVW in which the strain /uniFB01eld, denoted by b", is experimentally determined using full-/uniFB01eld 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) :b", the governing equation becomes Z [C(1;;k) :b"] :"?dV=Z SfTu?dS8u?2C(0). (/eight.pnum) Then, at least kindependent virtual /uniFB01elds must be chosen, among an in/uniFB01nite number of possibilities, to devise a system of kequations whose solution yields the kunknown constitutive parameters. The choice of these virtual /uniFB01elds is a key step. Three main approaches have been studied so far [Avril et al., /two.pnum/zero.pnum/zero.pnum/eight.pnum]: /one.pnum. manually using a polynomial function (see the example below); /two.pnum. automatically with special /uniFB01elds [Grédiac et al., /two.pnum/zero.pnum/zero.pnum/two.pnuma] that directly supply the constitutive param- eters from the virtual work of the external loads. These special /uniFB01elds also reduce the technique’s noise sensitivity [Avril et al., /two.pnum/zero.pnum/zero.pnum/four.pnum]; /three.pnum. piece-wise within the material [Toussaint et al., /two.pnum/zero.pnum/zero.pnum/six.pnum]. Continuous lower-degree polynomial func- tions between each sub-region can be used. These latter avoid to use a high-degree polynomial de/uniFB01ned on the whole geometry, and which would magnify the noise negative e/uniFB00ects. It should be noted that the volume integral in Eq. (/eight.pnum) requires the knowledge of the strain /uniFB01eld in the solid. This is a drawback of the VFM since most of full-/uniFB01eld 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 /uniFB01elds [Grédiac et al., /one.pnum/nine.pnum/nine.pnum/nine.pnum, Grédiac et al., /two.pnum/zero.pnum/zero.pnum/two.pnumb, Toussaint et al., /two.pnum/zero.pnum/zero.pnum/six.pnum]; /six.pnum •the properties of a non-linearly elastic material assuming a polynomial stress/strain relationship [Grédiac et al., /two.pnum/zero.pnum/zero.pnum/two.pnumb]; •the bending rigidities of anisotropic thin plates within the framework of the Kirchho/uniFB00-Love theory [Grédiac, /one.pnum/nine.pnum/eight.pnum/nine.pnum, Grédiac et al., /two.pnum/zero.pnum/zero.pnum/three.pnum]; •the sti/uniFB00ness and damping properties of thin viscoelastic isotropic vibrating plates [Giraudeau and Pierron, /two.pnum/zero.pnum/zero.pnum/five.pnum]; •the properties of a material within a plastic behavior using an iterative procedure based on a min- imization process as there is no direct relation between and the measured strains b"[Toussaint et al., /two.pnum/zero.pnum/zero.pnum/six.pnum]; •the damage parameters [Pierron and Grédiac, /two.pnum/zero.pnum/one.pnum/two.pnum] within the continuum damage mechanics [Lade- veze and LeDantec, /one.pnum/nine.pnum/nine.pnum/two.pnum]. 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 /uniFB01elds. /two.pnum./two.pnum Peridynamic (PD) framework /two.pnum./two.pnum./one.pnum Fundamental equations In peridynamics [Silling, /two.pnum/zero.pnum/zero.pnum/zero.pnum, Silling and Lehoucq, /two.pnum/zero.pnum/one.pnum/zero.pnum, Bobaru et al., /two.pnum/zero.pnum/one.pnum/six.pnum], a material point located atxinteracts with all other material points within a region Hx, centered on xand of a radius , which is referred to as the horizon (see Fig. /one.pnum). These material points in interaction with xare de/uniFB01ned in the reference con/uniFB01guration 0and are not updated over time. Let a material point be located at x0inHx. Letfv=fv(x0;x)be the dual force density that a point at x0 exerts on a point at x.fvis a force per unit volume squared and has the following properties (8x0=2Hx=)fv=0 (/nine.pnuma) fv(x;x0) =fv(x0;x) (/nine.pnumb) Eq. (/nine.pnumb)’s anti-symmetry results from the reciprocity principle. Thus, fvcan be expressed using the bond force densitytvas fv(x0;x) =tv(x0;x)tv(x;x0): (/one.pnum/zero.pnum) When assuming that long-range forces ( e.g., gravity) are negligible, and similarly to Eqs. (/one.pnum) and (/two.pnum), the PD e2Hx e1e3H x0 x 0 @ 0 Figure /one.pnum: A material point located at xinteracts with all other material points within a region H=Hx\ 0, centered onxand of a radius , which is referred to as the horizon [Delorme, /two.pnum/zero.pnum/one.pnum/eight.pnume] /seven.pnum equilibrium equation is 8 >>>>>>< >>>>>>:Z Hx\ 0fv(x0;x) dVx0+b(x) =08x2 (/one.pnum/one.pnuma) b=FonVf (/one.pnum/one.pnumb) u=uonVu (/one.pnum/one.pnumc) Fanduare the prescribed non-local boundary conditions acting on the boundary volume regions Vf andVu. These are thin real material layers on the boundary of and are such that Vf[Vu=@ and Vf\Vu=;.bis the external load applied as a body force density, i.e.,b=0in rVf. This is di/uniFB00erent from classical continuum mechanics in which the boundary conditions are prescribed over the surface @ . Fig. /one.pnum 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, /two.pnum/zero.pnum/zero.pnum/eight.pnum, Silling and Lehoucq, /two.pnum/zero.pnum/one.pnum/zero.pnum]. /two.pnum./two.pnum./two.pnum State-based theory A mathematical object called states was developed [Silling et al., /two.pnum/zero.pnum/zero.pnum/seven.pnum, Silling and Lehoucq, /two.pnum/zero.pnum/one.pnum/zero.pnum] for the state-based theory. The states and their properties are presented in details in App. A of [Delorme et al., /two.pnum/zero.pnum/one.pnum/seven.pnum]. The main de/uniFB01nitions are recalled here. A vector state is a general mathematical object that maps vectors onto vectors. The well-known 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 Tlinked to a deformation state Y. Letx2 0andx02Hxdenote two material points in the reference con/uniFB01guration. Let be the vector de/uniFB01ned by=x0xand called the bond connected tox. The familyHofxis de/uniFB01ned by H= 2 R3n0 j(+x=x0)2(Hx\ 0) (/one.pnum/two.pnum) whereHandHxare schematized in Fig. /one.pnum. LetYbe the deformation state that transforms bonds connected to xto their deformed images as 82H;Y[x]hi=y(x+)y(x) =y(x0)y(x)2R3(/one.pnum/three.pnum) wherey(x)andy(x0)are, respectively, the deformed positions of xandx0. Fig. /two.pnum illustrates these vari- ables. It is further assumed that 8(x1;x2)2 2 0jx16=x2=)y(x1)6=y(x2) =)Y[x]hi6=0; (/one.pnum/four.pnum) which means that two distinct points in the reference con/uniFB01guration 0are also distinct in the deformed con/uniFB01guration t. Thus, the deformed direction vector state Mis de/uniFB01ned by M=Y jYj. (/one.pnum/five.pnum) The force state Tis related to the bond force density tvthrough 82H;T[x]hi=T[x]hx0xi=tv(x0;x)2R3. (/one.pnum/six.pnum) Combining Eqs. (/nine.pnum), (/one.pnum/one.pnuma) and (/one.pnum/six.pnum) leads to the PD equilibrium equation such as, 8x2 0, 8 >>>< >>>:Z 0 T[x]hx0xiT[x0]hxx0i dVx0+b(x) =0 (/one.pnum/seven.pnuma) fv(x0;x) =T[x]hx0xiT[x0]hxx0i (/one.pnum/seven.pnumb) App. B of [Delorme et al., /two.pnum/zero.pnum/one.pnum/seven.pnum] presents the three PD constitutive models which are the bond-based , the ordinary state-based andnon ordinary state-based theories. In the sequel, only the ordinary state based approach is used since isotropic linearly elastic materials are considered. /eight.pnum /two.pnum./two.pnum./three.pnum /two.pnumD 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). /two.pnumD plane stress constitutive model 8 >>>>>>>>>>< >>>>>>>>>>:T=1 3 221 1 s++ 1 1 d wes+ dwed M (/one.pnum/eight.pnuma) s=1 q" 9K+G+ 1 212# (/one.pnum/eight.pnumb) d=1 q(8G) (/one.pnum/eight.pnumc) where sand dare the spherical and deviatoric PD parameters, esandedthe spherical and deviatoric extension state, the Poisson’s ratio and qthe PD weighted volume (see Eq. (/three.pnum/seven.pnum) and [Silling et al., /two.pnum/zero.pnum/zero.pnum/seven.pnum] for more details). /two.pnumD plane strain constitutive model 8 >>>>>>>>< >>>>>>>>:T=1 3 2 s d wes+ dwed M (/one.pnum/nine.pnuma) s=1 q(9K+G) (/one.pnum/nine.pnumb) d=1 q(8G) (/one.pnum/nine.pnumc) /two.pnum./two.pnum./four.pnum Discretization Several approaches have been used to discretize the PD equilibrium equation (/one.pnum/seven.pnum) such as Galerkin /uniFB01nite el- ement methods [Chen and Gunzburger, /two.pnum/zero.pnum/one.pnum/one.pnum], Gauss quadrature [Weckner and Emmrich, /two.pnum/zero.pnum/zero.pnum/five.pnum], the spatial discretization [Emmrich and Weckner, /two.pnum/zero.pnum/zero.pnum/seven.pnum, Parks et al., /two.pnum/zero.pnum/zero.pnum/eight.pnum] or the EMU (name of the /uniFB01rst PD software implemented [Littlewood, /two.pnum/zero.pnum/one.pnum/five.pnum]) nodal discretization (ND). The last method (EMU-ND) consists of a mid- point (or one-point) quadrature in a Lagrangian spatial discretization. The EMU-ND has been chosen for its e/uniFB03cient load distribution scheme and acceptable computation time. Fig. /three.pnum shows the mesh of a reference position in /two.pnumD with regularly spaced nodes where his the /uniFB01xed mesh size. Each node i, located atxiin the 0 tx x0 y(x) y(x0)Y[x]hiReference position Deformed positionHHY[x]hi Figure /two.pnum: The deformation state Ytransforms bonds 2H to their deformed images. y(x)andy(x0) are, respectively, the deformed positions of xandx0[Delorme, /two.pnum/zero.pnum/one.pnum/eight.pnumc]. /nine.pnum dual mesh, is associated to a /uniFB01xed volume Vi. Assume that the volumes surrounding the nodes are non overlapping ( Vi\Vj=;), and are recovering the reference volume V , meaning thatPn iVi=V . The discrete familyHiof a nodeiis de/uniFB01ned by: Hi=fjjxj2 0andkik2=kxjxik2g (/two.pnum/zero.pnum) Thus, the discrete PD equilibrium equation under the EMU-ND scheme, and which has to numerically be solved for each node i, reads 8 >>< >>:X j2H i T[xi]hxjxiiT[xj]hxixji Vj+bi=0 (/two.pnum/one.pnuma) fv(xj;xi) =T[xi]hxjxiiT[xj]hxixji (/two.pnum/one.pnumb) wherebiis the external load applied to the node i. Recall thatbi=0if nodeiis 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 jwhoseVjis partially inside the horizon of a node iand partially outside (see Fig. /three.pnum). Several algorithms were evaluated [Seleson, /two.pnum/zero.pnum/one.pnum/four.pnum] 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., /two.pnum/zero.pnum/zero.pnum/eight.pnum] which has been implemented in the code [Delorme et al., /two.pnum/zero.pnum/one.pnum/eight.pnum] presented in Sec. /five.pnum. /three.pnum The Virtual Fields Method within Peridynamics /three.pnum./one.pnum Principle of Virtual Works Since Eq. (/one.pnum/seven.pnum) is valid for any x2 0, each term can be “multiplied” (using the Euclidean inner product) by any continuous vector function u?such as 0 @Z 0fv(x0;x) dVx0+b(x)1 Au?(x) = 0 . (/two.pnum/two.pnum) Integrating Eq. (/two.pnum/two.pnum) over the whole domain 0yields the PVW stated as Z 0Z 0fv(x0;x)u?(x) dVx0dVx+Z 0b(x)u?(x) dVx= 0. (/two.pnum/three.pnum) 0 h e1e2xiHi Vi Figure /three.pnum: /two.pnumD EMU-ND and its mesh with a constant mesh size h.Viis the surrounding volume of the node ilocated atxiandHiits discrete family. The integration over Hxfrom Eq. (/one.pnum/one.pnum) is now replaced by the discrete sum fron Eq. (/two.pnum/one.pnum) over Hi[Delorme, /two.pnum/zero.pnum/one.pnum/eight.pnumd]. /one.pnum/zero.pnum Assume that the virtual quantity u?is kinematically admissible, i.e.u?2C(0), then, the contribution of the unknown force over the boundary volume constraint region Vuvanishes. Eq. (/two.pnum/three.pnum) can be divided into the following “physical” quantities: 8 >>>>>>< >>>>>>:W? int=Z 0Z 0fv(x0;x)u?(x) dVx0dVx (/two.pnum/four.pnuma) W? ext=Z 0b(x)u?(x) dVx (/two.pnum/four.pnumb) whereW? intis the virtual work done by the internal dual force densities and W? extthe virtual work done by the external forces. As for CCM version of the PVW, W? int+W? ext= 0. /three.pnum./two.pnum Virtual Fields Method and material identi/uniFB01cation The PD ordinary state-based constitutive equations for linearly isotropic elasticity depend non-linearly on the material properties (Eqs. (/one.pnum/eight.pnum) and (/one.pnum/nine.pnum)). Because of this non-linear relationship, the ksought material parameters are solution of a non-linear system. Rather than explicitly solving this non-linear system, the non-linear set of at least kindependent equations were solved by minimizing a residual function . Since linearly isotropic elastic materials are considered, requires at least two independent virtual /uniFB01elds u?(k) to properly identify the material parameters ( K;G ). Thus, the residual function , de/uniFB01ned as a normalized squared residual between W? intandW? ext, was (K;G) =vuuuuuuuutM2X k=1 W?(k) int+W?(k) ext2 M2X k=1 W?(k) ext2; (/two.pnum/five.pnum) Note that the normalization was used only for comparison purpose. Also, it was assumed that u?(k)are such thatX k W?(k) ext2 6= 0. /four.pnum Case study: three-point bending test /four.pnum./one.pnum Geometry and boundary conditions Fig. /four.pnum describes a virtual three-point bending setup in which a beam, of dimensions LWt(length widththickness) and simply supported on two supports spaced of S. The coordinate system ( O;e1;e2) 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 F(0;W) =Fe2 (/two.pnum/six.pnum) /four.pnum./two.pnum The virtual displacement /uniFB01elds de/uniFB01nition Two trivial independent continuous virtual displacement /uniFB01elds u?(1)andu?(2)were de/uniFB01ned: u?(1)(x1;x2) = jx1jS 2 e2 (/two.pnum/seven.pnum) u?(2)(x1;x2) = jx1jx1x2 LW e1 (/two.pnum/eight.pnum) These virtual /uniFB01elds, depicted in Fig. /five.pnum, are kinematically admissible, i.e.they vanish on the supports: u?(1) S 2;0 =u?(2) S 2;0 =0. (/two.pnum/nine.pnum) /one.pnum/one.pnum /four.pnum./three.pnum Application of the Virtual Field Method The internal virtual works W? intwere approximated from Eqs. (/two.pnum/one.pnum) and (/two.pnum/four.pnuma) as W? int=X iX j2H i fv(xj;xi)u? i VjVi. (/three.pnum/zero.pnum) The external virtual works W? extwere calculated from Eq. (/two.pnum/four.pnumb) and are equal to: 8 >>< >>:W?(1) ext=Fu?(1) A=Fu?(1)(0;W) =FS 2(/three.pnum/one.pnuma) W?(2) ext=Fu?(2) A=Fu?(2)(0;W) = 0 (/three.pnum/one.pnumb) Finally,KandGwere extracted using the residual function from Eq. (/two.pnum/five.pnum) and the minimization process described in Sec. /five.pnum. F e1e2 LSWA Figure /four.pnum: Three-point bending setup in which a beam, of dimensions LWtand simply supported on two supports spaced of S, is submitted to an external force F. The coordinate system ( O;e1;e2) origin is located at the specimen’s bottom face center [Delorme, /two.pnum/zero.pnum/one.pnum/eight.pnuma]. F e1e2 (a) Shape of the /uniFB01rst virtual /uniFB01eld u?(1) F e1e2 (b) Shape of the second virtual /uniFB01eld u?(2) Figure /five.pnum: Virtual displacements applied to the three-point bending test. The coordinate system ( O;e1;e2) origin is located at the specimen’s bottom face center [Delorme, /two.pnum/zero.pnum/one.pnum/eight.pnumf]. /one.pnum/two.pnum /five.pnum Minimization process implementation /five.pnum./one.pnum Algorithm Fig. /six.pnum shows the algorithm state chart diagram used for the minimization process to extract the material properties from the actual displacement /uniFB01eld. The algorithm has in input the initial guess for the material properties (K;G) and the actual displacement /uniFB01eld bu. A uniformly distributed random value in a search intervalIsu/uniFB03ciently large to the include material’s expected properties was chosen as the initial guess. The gathered steps in the box (see Fig. /six.pnum) are repeated until the residual is smaller or equal to a /uniFB01xed tolerance. The main steps are described below: /one.pnum. The dual force density fvis computed from Eq. (/two.pnum/one.pnum) using the trial material properties and bu. For this step, the in-house state-based PD code PeriPyDIC [Delorme et al., /two.pnum/zero.pnum/one.pnum/eight.pnum] is used with the /two.pnumD constitutive modeling described in Sec. /two.pnum./two.pnum./three.pnum). Note that instead of solving a direct problem in which the displacement /uniFB01eld is computed from the boundary conditions, here the dual force density fvis evaluated from a provided displacement /uniFB01eld. /two.pnum. The internal virtual works W? intare computed from Eq. (/three.pnum/zero.pnum) with the in-house code PeriPyVFM [De- lorme and Diehl, /two.pnum/zero.pnum/one.pnum/eight.pnum], with the two virtual /uniFB01elds given in Eqs. (/two.pnum/seven.pnum) and (/two.pnum/eight.pnum). Note that PeriPyVFM is based on the discretized version of the equations given in Sec. /three.pnum and can be extended to as many virtual /uniFB01elds as required. /three.pnum. The residual de/uniFB01ned in Eq. (/two.pnum/five.pnum) is computed within PeriPyVFM usingW? intand the external virtual worksW? extcalculated from Eq. (/three.pnum/one.pnum). As long as the residual is greater than the tolerance , the “black box” optimization solver NOMAD [Le Di- gabel, /two.pnum/zero.pnum/one.pnum/one.pnum], 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 /uniFB01nds suitable values for the material properties and the algorithm terminates. Table /one.pnum summarizes the versions of the software used in this paper. Start Initial guess (K;)2Rand actual displacement /uniFB01eld bu PeriPyDIC : Compute the dual force density fv(Eq. (/one.pnum/seven.pnum)) PeriVFM : Compute the internal works W? int(Eq. (/three.pnum/zero.pnum)) PeriVFM : Compute the residual (Eq. (/two.pnum/five.pnum))  NOMAD : Get new (K;) Material properties (K;)foundYesNo Figure /six.pnum: State chart diagram of the algorithm to extract the material properties from the actual displace- ment /uniFB01eld. The non-linear black box optimization software NOMAD [Le Digabel, /two.pnum/zero.pnum/one.pnum/one.pnum], based on the MADS algorithm, is used for the minimization process. /one.pnum/three.pnum /five.pnum./two.pnum Choice of the horizon  Choosing the ratio =hof horizon and nodal spacing to gain convergence is still an open question within the PD community [Bobaru et al., /two.pnum/zero.pnum/zero.pnum/nine.pnum, Bobaru and Hu, /two.pnum/zero.pnum/one.pnum/two.pnum]. Local limits and asymptotically compatible discretizations were considered in [Du and Tian, /two.pnum/zero.pnum/one.pnum/five.pnum] and [Du, /two.pnum/zero.pnum/one.pnum/six.pnum]. It is shown in [Tian and Du, /two.pnum/zero.pnum/one.pnum/three.pnum] that, if the nodal spacing h!0decays faster than the horizon !0, the EMU-ND converges to the continuum local partial di/uniFB00erential equation (PDE) limit. An empirical choice of the horizon for a speci/uniFB01c material is to adjust the horizon to /uniFB01t experimental results. Convergence analysis can also be performed as described in [Bobaru et al., /two.pnum/zero.pnum/zero.pnum/nine.pnum] 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,mis the only parameter to adjust with respect to the residual function . Thus, a three- point bending test simulated by /uniFB01nite element with a mesh-size of h= 0:25mmwas used to calibrate m. The residual was computed for the exact material properties for di/uniFB00erent value of m. Fig. /seven.pnum plots as a function ofm. It can be seen that is large form2, gets smaller for the next values and is the smallest form= 7. Then, increases for the next two values. This kind of non stabilized behavior has been shown in [Diehl et al., /two.pnum/zero.pnum/one.pnum/six.pnum]. Therefore, m= 7 was used in the sequel to extract the material properties ( K;G ) using the algorithm presented in Fig. /six.pnum. Note that this is not a convergence study but only a simplistic way to choose the horizon . Table /one.pnum: Software used for the implementation of the state-based PD-VFM. Name Version Reference NOMAD 3:8 [Le Digabel, /two.pnum/zero.pnum/one.pnum/one.pnum] PeriPyDIC 0:1 [Delorme et al., /two.pnum/zero.pnum/one.pnum/eight.pnum] PeriPyVFM 0:1 [Delorme and Diehl, /two.pnum/zero.pnum/one.pnum/eight.pnum] 1 2 3 4 5 6 7 8 9 10 m0.00.20.40.60.81.01.2Residual Φ Figure /seven.pnum: Residual for the di/uniFB00erent m-values with h= 0:25mm. The objective of this simplistic analysis is to /uniFB01nd the best m-value to use in the PD model in order to obtain the smallest residual for the prescribed mesh-size. /one.pnum/four.pnum /six.pnum Numerical validation /six.pnum./one.pnum Actual displacements simulated by Finite Element Analysis The virtual three-point bending test presented in Fig. /four.pnum was modelled in /uniFB01nite element analysis (FEA) with ANSYS®. The setup size parameters were: L= 128 mm,W= 32 mm,t= 12:7mm andS= 96 mm. The geometry was discretized using the PLANE/one.pnum/eight.pnum/two.pnum /two.pnum-D Structural Solid element with a mapped meshing of h= 0:25mm. The element was used as a plane stress element and was de/uniFB01ned by four nodes having two translations at each node in e1ande2directions. A null displacement condition was applied to the nodes located on the two supports and the beam was loaded with a force of magnitude Fof1 000 Napplied to nodeAshown in Fig. /four.pnum. The Young’s modulus Ewas of 4 000 MPa and the Poisson’s ratio of0:3, which yields a bulk modulus K= 3 333:33MPa and a shear modulus G= 1 538:46MPa . Figs. /eight.pnuma and /eight.pnumb respectively show the components u1andu2for the simulated displacement /uniFB01eld buFEA. /six.pnum./two.pnum Extracted material properties using the displacement /uniFB01eld The displacement /uniFB01eld buFEAwas input in the algorithm described in Sec. /five.pnum to identify the material prop- erties. The search interval for the sought (K;G )wasI= [1 000; 10 000]2. This interval includes the entered values K= 3 333:33MPa andG= 1 538:46MPa . Based on Eq. (/three.pnum/one.pnum), the external virtual works wereW?(1) ext = 48 000 N mm andW?(2) ext= 0N mm . The extracted material properties are presented in Table /two.pnum. (a)u1-component of the displacement /uniFB01eld buFEA. (b)u2-component of the displacement /uniFB01eld buFEA. Figure /eight.pnum: Components u1(a) andu2(b) of the simulated displacement /uniFB01eld buFEA. /one.pnum/five.pnum When the initial guesses for KandGwere 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 /one.pnum/four.pnum more times with di/uniFB00erent initial guesses distributed within the search interval to investigate the method’s robustness and stability. Table /two.pnum gathers the results and shows that all the extracted material properties converge to the same values (they only di/uniFB00er for the fourth decimal point), even though the number of evaluations was di/uniFB00erent. Therefore, the minimization of the residual seems to be stable and means that exactly one global minimum must exist within the search interval. /six.pnum./three.pnum Noise analysis The actual displacement /uniFB01eld buFEAused to identify the properties was optimal with respect to bias and noise as this is a simulated /uniFB01eld. Thus, the slight di/uniFB00erences 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 e/uniFB00ect [Le and Bobaru, /two.pnum/zero.pnum/one.pnum/seven.pnum], the numerical round o/uniFB00 errors or the virtual /uniFB01eld chosen which can be more or less sensitive to the previous errors e/uniFB00ects [Pierron and Grédiac, /two.pnum/zero.pnum/one.pnum/two.pnum]. Real displacement /uniFB01elds are disturbed both by bias and noise [Sutton et al., /two.pnum/zero.pnum/zero.pnum/nine.pnum]. Contrary to the DIC bias that can be reduced or eliminated with a proper setup and parameters (good calibration, no contami- nation 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 buFEAas buNOISE=buFEA+gNOISE (/three.pnum/two.pnuma) gNOISE=g1e1+g2e2 (/three.pnum/two.pnumb) wheregNOISE is the disturbing displacement vector whose components g1andg2follow a Gaussian distri- bution of mean and standard deviation . It was assumed that = 0:0mm, and lied between 0:0 and2:5104mmwhich is a typical standard deviation range measured with the DIC-system used in the Laboratory for Multiscale Mechanics (LM/two.pnum). Table /two.pnum: Extracted material properties from buFEA.15initial guesses (K;G)2[1 000; 10 000]2were used and provided in the second and third columns. The fourth and /uniFB01fth columns represent the extracted ma- terial 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 #K(MPa )G(MPa )K(MPa )G(MPa )K(%)G(%) Eval /one.pnum5 000:000 5 000 :000 3 354:194 1 531 :635 0:630:44 /four.pnum/three.pnum/five.pnum /two.pnum1 098:722 7 031 :443 3 354:194 1 531 :635 0:630:44 /five.pnum/zero.pnum/eight.pnum /three.pnum2 495:136 6 502 :603 3 354:194 1 531 :635 0:630:44 /four.pnum/nine.pnum/three.pnum /four.pnum9 845:479 1 797 :806 3 354:194 1 531 :635 0:630:44 /five.pnum/zero.pnum/nine.pnum /five.pnum1 858:444 9 782 :045 3 354:194 1 531 :635 0:630:44 /five.pnum/three.pnum/two.pnum /six.pnum4 334:080 5 263 :272 3 354:194 1 531 :635 0:630:44 /five.pnum/one.pnum/zero.pnum /seven.pnum6 410:611 3 680 :783 3 354:194 1 531 :635 0:630:44 /four.pnum/five.pnum/zero.pnum /eight.pnum9 881:218 2 116 :642 3 354:194 1 531 :635 0:630:44 /three.pnum/nine.pnum/zero.pnum /nine.pnum8 723:166 7 719 :498 3 354:194 1 531 :635 0:630:44 /four.pnum/two.pnum/seven.pnum /one.pnum/zero.pnum 3 718:518 6 977 :810 3 354:194 1 531 :635 0:630:44 /four.pnum/three.pnum/one.pnum /one.pnum/one.pnum8 486:055 9 008 :744 3 354:194 1 531 :635 0:630:44 /five.pnum/four.pnum/one.pnum /one.pnum/two.pnum 3 396:443 6 769 :252 3 354:194 1 531 :635 0:630:44 /five.pnum/three.pnum/eight.pnum /one.pnum/three.pnum 5 543:759 3 163 :586 3 354:194 1 531 :635 0:630:44 /five.pnum/one.pnum/five.pnum /one.pnum/four.pnum 5 723:661 6 488 :884 3 354:194 1 531 :635 0:630:44 /four.pnum/three.pnum/five.pnum /one.pnum/five.pnum 4 612:369 1 874 :837 3 354:194 1 531 :635 0:630:44 /five.pnum/four.pnum/three.pnum /one.pnum/six.pnum Table /three.pnum provides the extracted material properties after adding Gaussian noise. The relative error remained under 1:4%in absolute value. Therefore, the virtual /uniFB01elds chosen seem stable with respect to noise. This could be di/uniFB00erent with other virtual /uniFB01elds. Also, it should be noted that the bulk modulus Kis more sensitive to Gaussian noise than the shear modulus G. The conclusion from this noise analysis is similar that presented in [Kramer and Scherzinger, /two.pnum/zero.pnum/one.pnum/four.pnum], in which the VFM is applied within CCM. /six.pnum./four.pnum Missing data Obtaining displacement /uniFB01elds on the exact specimen edges is hardly possible with DIC and the missing data due to this experimental limitation could have a signi/uniFB01cant in/uniFB02uence on the identi/uniFB01cation results, depending the the virtual /uniFB01elds chosen [Pierron and Grédiac, /two.pnum/zero.pnum/one.pnum/two.pnum]. Lines of elements that were placed of a multiple of h= 0:25mmhave been removed from the specimen’s edges (see Fig. /nine.pnum) before extracting the material properties to simulate the e/uniFB00ect of the missing data. Fig. /one.pnum/zero.pnum presents the e/uniFB00ects 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 /four.pnum provides the error in the identi/uniFB01cation of KandG, showing that the relative errors quickly increase for the bulk modulus K. This means the identi/uniFB01cation of Kis more sensitive to missing data than the shear modulus Gusing the three-point bending experimental setup and the virtual /uniFB01elds de/uniFB01ned in Eqs. (/two.pnum/seven.pnum) and (/two.pnum/eight.pnum). /seven.pnum Conclusion The main contributions of this study are as follows: /one.pnum. The VFM within PD has been derived to extract the elastic properties of a linearly isotropic material. Open-source python codes [Delorme and Diehl, /two.pnum/zero.pnum/one.pnum/eight.pnum,Delorme et al., /two.pnum/zero.pnum/one.pnum/eight.pnum] 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 bene/uniFB01ts of the PD-VFM is that no computation of the displacement /uniFB01eld during the minimization process is required. /two.pnum. The PD-VFM has been tested on a three-point bending test using undisturbed displacement /uniFB01elds Table /three.pnum: Error in the identi/uniFB01cation of KandGdue to Gaussian noise. The /uniFB01rst column is the noise’s stan- dard deviation used. The second and third columns represent the extracted material properties after minimization. The fourth and /uniFB01fth 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) Standard Deviation Minimization results Relative error (mm)K(MPa )G(MPa )K(%)G(%) 0:01043 354:194 1 531:635 0:630:44 0:51043 347:890 1 531:918 0:440:43 1:01043 379:440 1 531:790 1:380:43 1:51043 339:325 1 531:485 0:180:45 2:01043 323:665 1 532:6920:290:37 2:51043 364:742 1 533:479 0:940:32 Missing nodeh Figure /nine.pnum: Missing nodes all around the specimen’s edges [Delorme, /two.pnum/zero.pnum/one.pnum/eight.pnumb]. /one.pnum/seven.pnum generated by FEA and two independent virtual /uniFB01elds 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 prop- erties identi/uniFB01cation within less than /zero.pnum./seven.pnum/five.pnum% 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 mea- surements performed on a polycarbonate beam under a three-point bending test. However, the residual computed using the actual material properties (measured per ASTM D/six.pnum/three.pnum/eight.pnum-/one.pnum/four.pnum) was too high and then did not allow to extract the material properties. This can be explained by the di/uniFB03culty of reducing the missing data due to the DIC measurements and the high sensibility of the bulk modulus Kwith respect to this miss- ing information and virtual /uniFB01elds chosen. Based on the full-/uniFB01eld measurement literature [Belhabib et al., /two.pnum/zero.pnum/zero.pnum/eight.pnum,Pierron and Grédiac, /two.pnum/zero.pnum/one.pnum/two.pnum], the authors would recommend to design a well-suited test con/uniFB01guration for the PD-VFM that would meets the following criteria: 0.00 0.25 0.50 0.75 1.00 1.25 Missing distance (mm)0.000.010.020.030.040.05Residual Φ Figure /one.pnum/zero.pnum: In/uniFB02uence of missing data: the black line shows the linear regression between the residual with respect to the missing distance around the specimen. Table /four.pnum: Error in the identi/uniFB01cation of KandGdue to the missing data on the specimen free edges. The /uniFB01rst column is the missing distance on the edge (see Fig. /nine.pnum). The second and third columns represent the extracted material properties after minimization. The fourth and /uniFB01fth 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) Missing distance Minimization results Relative error mm K(MPa )G(MPa )K(%)G(%) 0:00 3 354:194 1 531:635 0:630:44 0:25 3 577:209 1 552:676 7:32 0:92 0:50 3 708:762 1 565:134 11:26 1:73 0:75 3 853:462 1 578:672 15:60 2:61 1:00 4 011:848 1 593:237 20:36 3:56 1:25 4 184:931 1 608:811 25:55 4:57 /one.pnum/eight.pnum •a large heterogeneity of the displacement /uniFB01elds. This is the spirit of the full-/uniFB01eld measurement meth- ods, by opposition to standard mechanical tests that are based on homogeneous and averaged /uniFB01elds; •a good sensitivity of the displacement /uniFB01elds with respect to the material properties to identify. The authors also point out that they were only interested in this study to present a proof of concept of the PD-VFM by extracting the elastic properties of a linearly isotropic material. However, they believe that fu- ture works could use the approach presented through this paper and extend it for damage characterization, which is the main strength of PD. A Isotropic elastic materials within state-based peridynamics A./one.pnum Elastic material An elastic material is de/uniFB01ned by its free energy density which only depends on Y. De/uniFB01ne = (Y)as thestrain energy density . Then, the force state Tis de/uniFB01ned by T=T(Y) =@ @Y=d dY=O (Y). (/three.pnum/three.pnum) A./two.pnum Isotropy A material is isotropic if and only if, for any proper orthogonal second-order tensor Q T(YQ)hi=T(Y)hQi. (/three.pnum/four.pnum) Physically, Eq. (/three.pnum/four.pnum) re/uniFB02ects the fact that the force state is invariant by applying rotations before deforma- tions. A./three.pnum Isotropic elastic solid Let#be a non-local volume dilatation [Silling and Lehoucq, /two.pnum/zero.pnum/one.pnum/zero.pnum,Le et al., /two.pnum/zero.pnum/one.pnum/four.pnum,Sarego et al., /two.pnum/zero.pnum/one.pnum/six.pnum,Le and Bobaru, /two.pnum/zero.pnum/one.pnum/seven.pnum] such as #(e) =% q(wx)e (/three.pnum/five.pnum) wis the scalar in/uniFB02uence state used to weight the bonds in/uniFB02uence on the force state calculation and x= jXhij=kk2.qis the weighted volume such as q= (wx)x (/three.pnum/six.pnum) %=8 >>>< >>>:3 in /three.pnumD 221 1 in /two.pnumD plane stress 2 in /two.pnumD plane strain(/three.pnum/seven.pnum) is the Poisson’s ratio. Finally, eis the scalar extension state de/uniFB01ned by e=jYhijx=yx. (/three.pnum/eight.pnum) Thescalar extension state can be divided into a spherical part esand a deviatoric part edsuch as e=es+ed=x 3#(e) +ed. (/three.pnum/nine.pnum) Thespherical part represents the isotropic expansion while the deviatoric part depicts shear deforma- tions. /one.pnum/nine.pnum Suppose an isotropic elastic material de/uniFB01ned by its strain energy density (e) =1 2 s(wes)es+1 2 d(wed)ed(/four.pnum/zero.pnum) where sand dare the PD paremeters [Delorme et al., /two.pnum/zero.pnum/one.pnum/seven.pnum]. They are calibrated to the classical isotropic elastic material properties ( KandG, 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 in/uniFB01nitesimal strain theory [Silling et al., /two.pnum/zero.pnum/zero.pnum/seven.pnum]. One /uniFB01nds [Le et al., /two.pnum/zero.pnum/one.pnum/four.pnum,Sarego et al., /two.pnum/zero.pnum/one.pnum/six.pnum,Le and Bobaru, /two.pnum/zero.pnum/one.pnum/seven.pnum,Delorme et al., /two.pnum/zero.pnum/one.pnum/seven.pnum] that the PD parameters are: s=8 >>>>>>>< >>>>>>>:1 q(9K) in /three.pnumD 1 q" 9K+G+ 1 212# in /two.pnumD plane stress 1 q(9K+G) in /two.pnumD plane strain(/four.pnum/one.pnum) d=8 >>< >>:1 q(15G)in /three.pnumD 1 q(8G) in /two.pnumD plane stress or plane strain(/four.pnum/two.pnum) The calibration above also assumes that the material points have a complete neighborhood of size 2(i.e. H=Hx) /uniFB01lled with the same material. This assumption is not veri/uniFB01ed within a distance of next to a free surface or an interface with another materials and implies that Eqs. (/four.pnum/one.pnum) and (/four.pnum/two.pnum) are no more valid near boundaries. A corrective approach, called force normalization , was devised to account for that phe- nomenon [Macek and Silling, /two.pnum/zero.pnum/zero.pnum/seven.pnum]. From Eqs. (/three.pnum/eight.pnum), (/three.pnum/five.pnum), (/three.pnum/nine.pnum) and (/four.pnum/zero.pnum), one obtains: e= es+ ed(/four.pnum/three.pnuma) #(e) =% q(wx)e+o(jjejj) (/four.pnum/three.pnumb) es(e) =x 3% q(wx)e+o(jjejj) (/four.pnum/three.pnumc) ed(e) = 1x 3% q(wx) e+o(jjejj) (/four.pnum/three.pnumd)  (e) = s(wes)es+ d(wed)ed+o(jjejj) (/four.pnum/three.pnume) where 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 incremen- tal change to the deformation state Yas follows e=Y jYjY+o(kYk) =MY+o(kYk). (/four.pnum/four.pnum) Thus, combining Eqs (/three.pnum/eight.pnum), (/three.pnum/three.pnum), (/four.pnum/three.pnum) and (/four.pnum/four.pnum) leads to the constitutive model of an isotropic elastic material given by the force state T(Y) =1 3 % s(3%) d wes+ dwed M. (/four.pnum/five.pnum) Eq. (/four.pnum/five.pnum) shows that the constitutive model of an isotropic elastic material corresponds to an ordinary ma- terial since T=tM. /two.pnum/zero.pnum