Nonlocal strong forms of thin plate, gradient elasticity, magneto-electro-elasticity and phase field fracture by nonlocal operator method

The derivation of nonlocal strong forms for many physical problems remains cumbersome in traditional methods. In this paper, we apply the variational principle/weighted residual method based on nonlocal operator method for the derivation of nonlocal forms for elasticity, thin plate, gradient elasticity, electro-magneto-elasticity and phase field fracture method. The nonlocal governing equations are expressed as integral form on support and dual-support. The first example shows that the nonlocal elasticity has the same form as dual-horizon non-ordinary state-based peridynamics. The derivation is simple and general and it can convert efficiently many local physical models into their corresponding nonlocal forms. In addition, a criterion based on the instability of the nonlocal gradient is proposed for the fracture modelling in linear elasticity. Several numerical examples are presented to validate nonlocal elasticity and the nonlocal thin plate .


Introduction
Classical continuum mechanics has achieved great success in describing the macro-scale properties of solid material based on the continuous medium hypothesis that the material is a continuous mass rather than as discrete particles. The assumption indicates that the substance of the object completely fills the space it occupies, without considering the inherent micro-structure of the material. Such a continuous medium hypothesis is not always valid in solid medium. Over the years, researchers found that many phenomena, such as size effect [1], length scale effect [2], skin/edge effect [3], can not be well predicted by traditional continuum mechanics. These phenomena may be attributed to the nonlocal effect in the solid. In contrast with local theory whose mathematical language is partial differential derivatives defined at an infinitesimal point, nonlocal theory is formulated as integral form in a domain.
Classical continuum mechanics is regarded as a local theory. For solid mediums of multiple materials with material interface or discontinuity such as fracture, the partial differential operator is no longer well defined. Around the fracture front tip, the stress singularity happens for local theory. In order to model fracture and its evolution, various local theories have been proposed, for example, finite element method (FEM) [4], extended finite element method [5], phase-field fracture method [6,7,8], cracking particle method [9,10], extended finite element method [11], numerical manifold method [12], extended isogeometric analysis (XIGA) for three-dimensional crack [13], meshfree methods [14,15,16]. Another approach for fracture modeling is the nonlocal method. Compared with continuum mechanics without length scale, nonlocal theory takes into account the length scale explicitly and it is less sensitive to the inhomogeneity/discontinuity encountered in the materials due to its integral form.
Two general theories to account for the length scale of solid material, are the gradient elasticity [17,1,18,19] and the nonlocal elasticity [20,21,22,23]. The gradient elasticity theory can be traced back to Cosserat theory in 1909 [24]. It incorporates the length scale and higher order derivative of the displacement field. A variety of gradient elasticity theories have been proposed such as Mindlin solid theory [17,2], couple stress theory [1,25], modified couple stress [18,26] and second-grade materials [19]. In nonlocal elasticity, the stress tensor is based on the integral of the "local" stress field in a domain, in contrast with the local elasticity defining the stress based on the strain field at a point. Under certain circumstances, the nonlocal elasticity can be transformed into gradient elasticity [22,27].
Among various nonlocal elasticity theories, Peridynamics (PD) [28,29] have attracted the attention of the researchers in the fracture mechanic field. PD is based on the integral form well defined in domain with/without discontinuity. This salient feature enables PD a versatile method for fracture modeling [30,31,32,33]. The origin of PD is the bond-based PD (BB-PD) with the Poisson ratio restriction. BB-PD can model 2D elasticity with Poisson ratio of 1/3 and 3D elasticity with Poisson ratio of 1/4. Many efforts have been dedicated to overcome this restriction, for example, PD with shear deformation [34], bond-rotation effect by [35], PD with micropolar deformation [36]. The further development of PD is the state-based PD [29,37]. Several treatments are developed to overcome the instability issue in non-ordinate state-based PD (NOSBPD), including, bond-associated higher-order stabilized model [38], higher-order approximation [39], stabilized non-ordinary state-based PD [40,41], sub-horizon scheme [42] and stress point method [43].
Dual-horizon PD overcomes the restriction of constant horizon in PD, without introducing side effects for variable support size. Dual-horizon peridynamic formulation can be derived from the Euler-Lagrange equations [58]. Based on the concept in nonlocal theory, we developed the Nonlocal Operator Method (NOM) as the generalization of dual-horizon PD. NOM uses the nonlocal operators of integral form to replace the local partial differential operators of different orders. There are three versions of NOM, firstorder particle-based NOM [59,60], higher order particle-based NOM [61] and higher order NOM based on numerical integration [62]. The particle-based version can be viewed as a special case of NOM with numerical integration when nodal integration is employed. The nonlocal operators can be viewed as an alternative to the partial derivatives of shape functions in FEM. Combined with a variational principle or weighted residual method, NOM obtains the residual vector and tangent stiffness matrix in the same way as in FEM. NOM has been applied to the solutions of the Poisson equation in high dimensional space, von-Karman thin plate equations, fracture problems based on phase field [61], waveguide problem in electromagnetic field [60], gradient solid problem [62] and Cahn-Hilliard equation [63].
Although much progress in nonlocal methods has been achieved in the above mentioned literatures, the derivations for many physical problems remain cumbersome and complicated, see for example [47,64,57,65]. In local theory, the local differential operator is a fundamental element for describing physical problems. In analogy, the nonlocal operators would be very beneficial for developing nonlocal theoretical models. The power of NOM in deriving nonlocal models remains largely unexplored. In addition, NOM based on implicit algorithms is relatively complicated in implementation and in this paper, we explore the explicit algorithm in solving the nonlocal models. Furthermore, we propose an instability criterion of the nonlocal gradient operator for the purpose of fracture modeling. The remaining of the paper is outlined as follows. In section 2, the second-order NOM in 2D/3D is formulated in detail. In section 3, we apply the NOM scheme combined with variational principle/weighted residual method to derive the nonlocal governing equations for elasticity, thin plate, gradient elasticity, electro-magnetoelasticity and phase field fracture model. The correspondence between local form and nonlocal form for higher order problems is discussed. In section 4, an instability criterion of nonlocal gradient is presented in the fracture modeling of linear elastic solid. The implementation of nonlocal solid and nonlocal thin plate is discussed in section 5. Several numerical examples for solid and thin plate are used to demonstrate the accuracy and efficiency of the current method in section 6. Last but not the least, some concluding remarks are presented.
2 Second-order nonlocal operator method NOM uses the integral form to replace the partial differential derivatives of different orders. Although NOM can solve higher order linear/nonlinear problems in 2D/3D, we restrict our discussion in second-order NOM, which is sufficient for the nonlocal derivation of the physical problems to be studied in section 3.

Support and dual-support
Consider a domain as shown in Fig.1(a), let x i be spatial coordinates in the domain Ω; and v j := v(x j , t) are the field values for x i and x j , respectively; v ij := v j − v i is the relative field vector for spatial vector r. Support S i of point x i is the neighbourhood of point x i . A point x j in support S i forms the spatial vector r(= x j − x i ). The support in the NOM can be a spherical domain, a cube, semi-spherical domain and so on. Figure 1: (a) Domain and notation. (b) Schematic diagram for support and dual-support, all shapes above are supports, Dual-support is defined as a union of points whose supports include x i , denoted by Point x j forms the dual-vector r ji (= x i − x j = −r ij ) in S i . On the other hand, r ji is the spatial vector formed in S j . It is worth mentioning that the size of the support of each point can be different. When the support sizes for all material points are the same, the dual-support is equal to the support. On the other hand, if the size of support varies for each point, the shape of dual-support can be quite irregular, even discontinuous for two adjacent points. One example to illustrate the support and dual-support is shown in Fig.1(b).

Dual property of dual-support
For point j ∈ S i , let f ij be a physical quantity, work conjugate to field difference (u j − u i ), the dual property of dual-support is Proof : Let the domain Ω be divided into N non-overlapping particles, so that Ω = N i=1 ∆V i , where ∆V i is the volume assigned to particle i. Herein, N can be arbitrarily large so that the ∆V i is infinitesimal and the double summations of discrete form converge to the double integrals in continuous form.
In the third step, the dual-support is considered as follows. The term f ij with u j is the physical quantity from i's support, but is added to particle j; since j ∈ S i , i belongs to the dual-support S j of j; all terms f ji with u i are collected from any material point j whose support contains i and hence form the dual-support of i. Therefore, the dual property of the dual-support is proved.
When all points have the same size of support domains, i.e. j ∈ S i ↔ i ∈ S j , we have S i = S i for any point i and then the dual property of dual-support by Eq.2 becomes Above equation is widely used in the derivation of nonlocal strong form from weak form. Such expression is valid in the continuum form as well as in discrete form. The dual property of dual-support is also proved in the dual-horizon peridynamics [45]. A simple example with N = 4 to illustrate this property is given in A.

Nonlocal gradient and Hessian operator
The local gradient operator and Hessian operator for a scalar-valued function u have the forms in 2D ∇u = u ,x , u ,y T , ∇ 2 u = u ,xx u ,xy u ,xy u ,yy (5) and in 3D where u ,xx denotes the partial derivative of u with respect to x twice.
In the framework of NOM, the partial derivatives can be constructed as follows. The Taylor series expansion of scalar-valued field u j in 2D can be written as where r ij = (x ij , y ij ) T = x j − x i and O(|r ij | 3 ) denotes the higher order term. Let The Taylor series expansion of Eq.7 can be rewritten as Tensor product with p T ij on both sides of Eq.11 Considering the weighted integration in the support S i , we obtain where ω(r ij ) is the weight function.
Then the nonlocal operators can be obtained as where Here, we use˜ to denote the nonlocal form of the local operator since the definitions of the local operator and the nonlocal operator are distinct. The Taylor series expansion of a vector field u can be obtained in the similar manner as That is∂ For example, consider the displacement field u = (u, v) T in two dimensional space, the relative displacement vector and the nonlocal partial derivatives have the explicit forms Let K i · p ij be denoted by The gradient vector g ij and Hessian matrix h ij between points i and j in 2D are, respectively In 3D case, the polynomial vector based on relative coordinates The shape tensor in 3D is constructed by Eq.15 with p ij in Eq.23.
Let K i · p ij in 3D be denoted by The gradient vector g ij and Hessian matrix h ij for two points i, j in support in 3D are, respectively It is worth mentioning that for first order NOM or peridynamics, the gradient vector can be calculated as well by Then the nonlocal gradient operator and Hessian operator for vector field can be defined as∇ In the case of 2-vector in 2 dimensional space, the explicit forms of∇⊗u i and∇ ⊗∇ ⊗ u i are∇ For scalar-valued field, the nonlocal Laplace operator is the tensor contraction of∇ ⊗∇u i , e.g.∆ =∇ ·∇ = tr(∇ ⊗∇), where tr(·) denotes the trace of a matrix. More specifically, in 2D and in 3D∆ And their local counterparts for scalar-valued field are ∆w = w ,yy + 2w ,xy + w ,xx in 2D (33) ∆w = w ,xx + w ,yy + w ,zz + 2w ,xy + 2w ,xz + 2w ,yz in 3D (34)

Stability of the second-order nonlocal operators
According to Ref [61], the energy functional for second-order nonlocal operator in discrete form can be written as where p hg is the penalty and m i = S i ω(r ij ) dV j . The operator in Eq.14 corresponds to the minimum of Eq.35. The first variation of F i is For any δu i , Ω δF i dV i = 0 leads to the internal force due to the stability of the nonlocal operator Eq.38 is the expression for a scalar-valued field. For vector-valued field, the internal force due to the stability of nonlocal operator is

Nonlocal governing equations based on NOM
This section is devoted to the variational derivation of nonlocal strong forms of solid mechanics, including hyperelasticity, thin plate, gradient elasticity, electro-magnetic-elasticity theory and phase field fracture method. The strong form is suitable for theoretical analysis as well as explicit time integration. For the fully implicit simulation of various PDEs, the reader is referred to NOM for PDEs [59,60,61,62,63].

Nonlocal form for hyperelasticity
Consider the energy density of a hyperelasticity as φ := φ(F ), where F = ∇u + I. The balance equation for the hyperelastic solid is ∇ · P + b = 0 on Ω (40) with boundary conditions u = u 0 on Γ D and P · n = t 0 on Γ N , where u 0 is the specified displacement and t 0 is the prescribed traction load, P = ∂φ ∂F , the first Piola-Kirchhoff stress, b is the body force density.

Derivation based on variational principle
The variation of strain energy over the domain is In above derivation, we replace the gradient operator with nonlocal gradient, e.g.∇ ⊗ u i → S i ω(r ij )u ij ⊗ g ij dV j in Eq.27, and the relation A : a ⊗ b = (A · b) · a for second-order tensor A and vectors a, b is employed.
The variational of external body force energy For any δu i , δF − δF ext = 0 leads to the nonlocal governing equations for elasticity Considering the effect of inertial force ρü i per unit volume, and replacing the dual-support with dual-horizon, we obtain the equations of motion for dual-horizon peridynamics If the sizes of horizons for all material points are the same, the dual-horizon peridynamics degenerates to the conventional constant horizon peridynamics.
For any specific strain energy density (for example, isotropic/anisotropic linear/nonlinear elasticity), the explicit form of P can be derived straightforwardly.

Derivation based on weighted residual method
Beside the derivation based on strain energy density, the nonlocal strong form can be derived by weighted residual method. Consider the governing equations for hyperelasticity , the weak form of Eq.40 for any trial vector becomes Let us focus on the integral in Ω, the first term in above equation can be written as For any v i , the weak form being zero leads to which is identical to Eq.43. As being more general than the energy method, the weighted residual method can be used to convert PDEs that have no energy functional to nonlocal integral forms.

Nonlocal thin plate theory
The thin plate theory is widely used in engineering applications [66]. The basic assumption of thin plate include: 1) the thickness of the plate is much smaller than the length inside the mid-plane; 2) the deflection is much smaller than the thickness of the plate so that higher order effect is neglect-able; 3) the stress along the thickness direction is assumed as zero, e.g. σ z ≈ 0 and the points in the midplane have no displacement parallel to the midplane, e.g. u(x, y, 0) = v(x, y, 0) ≈ 0; 4) the normal of the mid-plane remains perpendicular to the mid-plane after deformation. Then the plate bending can be simplified into 2D problem and the displacements, strain and stress can be described by the deflection on the mid-plane The generalized strain is the Hessian operator on the deflection with nonlocal correspondence and its variation The momentum tensor M , the general stress for isotropic thin plate, is given by where D 0 = Et 3 12(1−ν 2 ) and t is the thickness of the plate. Based on the principle of minimum potential energy, the energy functional for the governing equation is (54) and for the boundary condition can be expressed as where q is the external transverse load on the mid-plane,V n is the shear force load on boundary S 3 andM n is the prescribed moment on boundary S 2 + S 3 . For simplicity, we leave the integral on the boundary for later consideration. The variation of the internal energy functional is The variation of the external energy function is For any δw i , δF int − δF ext = 0 leads to the nonlocal thin plate equation for material point in domain Ω The additional nonlocal form for material point applied with the moment boundary condition is Based on the D'Alembert's principle, the equation of motion considering the effect of inertial force ρtẅ i per unit area is For clamped boundary condition w ,n = ∇w · n = 0, the nonlocal form is Compared with the local governing equation for thin plate ∇ 2 : M + q = tρẅ, we can find the correspondence between local and nonlocal formulation The nonlocal derivation for thin plate can be extended to composite plate and functional gradient plate theories.

Nonlocal gradient elasticity
Gradient theories emerge from considerations of the microstructure in the material at micro-scale, where a mass point after homogenization is not the center of a micro-volume and the rotation of the micro-volume depends on the moment stress/couple stress as well as the Cauchy stress. Gradient elasticity generalizes the elasticity theory by employing higher order terms of the deformation gradient or the gradient of the strain tensor. Generally, the energy density functional can be assumed as ψ := ψ(F , ∇F ) = ψ(∇u, ∇ 2 u), where F = ∇u + I. The total potential energy in domain is The stress tensor and generalized stress tensor of first Piola-Kirchhoff type are defined as The variation of the total internal energy is δF = Ω ∂ψ ∂F : ∇δu + ∂ψ ∂∇F˙: Based on the integration by parts, the local form can be derived by Based on D'Alembert's principle, the governing equations for dynamic gradient elasticity can be written as On the other hand, do the substitutions ∇δu → S i ω(r ij )g ij ⊗δu ij dV j , and ∇ 2 δu → In the above derivation, we used Σ:u ⊗ h = (Σ : h) · u. For any δu i , δF = 0 leads to the nonlocal form of gradient elasticity The inertia force term is added based on D'Alembert's principle. Comparing Eq.67 and Eq.69, the correspondence from local form to nonlocal form is

Nonlocal form of magneto-electro-elasticity
In accordance with reference [67], let us postulate the following form of internal energy for the energy function ψ := ψ(F , ∇F , p, ∇p, m, ∇m), a function depends on the displacement gradient F = ∇u + I and its second gradient ∇F = ∇ 2 u, polarization vector p and its gradient ∇p, magnetic field m and its gradient ∇m. The total potential energy in the domain can be written as This model has a strong physical background, for example, the nonlinear electro-gradient elasticity for semiconductors [68] and flexoelectricity [69].
The first variation of F is δm ij ⊗g ij dV j and following the same operations in prior sections, the functional becomes For any δu i , δp i , δm i , δF = 0 leads to general nonlocal governing equation for mechanical field, electrical field and magnetic field, respectively In the derivation, we did not specify the exact form of the energy density, whether it is of small deformation or of finite deformation. For the specified energy form, one only needs to derive the expression for σ, Σ, e, E, s, S based on the material constitutions. It can be seen that the nonlocal governing equations for the continuum magneto-electro-elasticity can be obtained with ease by using nonlocal operator method and variational principle. The same rule applies for many other physical problems.

Nonlocal form of phase field fracture method
Phase field fracture method is powerful in fracture modelling [70]. The difference in tensile and compressive strengths of the material can be considered by dividing the strain energy density into a tensile part affected by the phase field and a compressive part, which is independent of the phase field, ψ e (ε(∇u), s) = (1 − s) 2 ψ + e (ε(∇u)) + ψ − e (ε(∇u)).
where ψ + e (ψ − e ) denotes the strain energy density for tensile (compressive) part, u the displacement, s ∈ [0, 1] the phase field, ε the strain and is the phase field intrinsic length scale.
The full potential functional of the phase field fracture model reads where t * the surface traction at the boundary, b the body force density and g c is the critical energy release rate.
For the sake of simplicity, we neglect the surface traction force and consider the first variation of F where For any δu i , δs i , δF = 0 leads to the nonlocal governing equations for the mechanical field and phase field The above examples aim at illustrating the power of nonlocal operator method combined with weighted residual method or variational principle in the derivation of nonlocal strong forms based on their local strong or energy forms. The derived nonlocal strong forms are variationally consistent and allow variable support sizes for each point in the model.

Instability criterion for fracture modelling
Typical methods for fracture modelling are either based on diffusive crack domain in phase field methods or on direct topological modification on meshes in XFEM or bonds in PD. Direct topological modification on meshes often leads to instability issues. For example, in NOSBPD, the breakage of a bond based on the quantities derived from stress state or strain state often introduces too much perturbation to the scheme, which may abort the calculation because of the singularity in shape tensors. These criteria include critical stretch [28,71], energy based [30] or stress based criterion [32,33]. Another issue in NOSBPD is that the strain energy carried by a bond is not independent with other bonds. It also depends on the direction, the length of the bond, the choice of influence functions. Removing one neighbour often gives rise to catastrophic result on the calculation. A criterion on how to remove the neighbours safely from the neighbour list remains unclear.
Damage is a process deviated from the robust mathematical expression, where the transition happens in a very narrow zone, such as the crack tip front. It is observed that around the crack tip, the gradient or strain undergoes a sharp transition within a very small zone. Most conventional numerical methods for fracture modelling focus on accurate description of the singularity occurring around the crack tip, such a description is very hard to tackle and its evolution is inconvenient to update. This dilemma can be handled when something different from continuous function is introduced.
In NOM, the gradient operator is defined in a "redundant" way. Around the crack tip, the deformation is irregular and the part due to hourglass energy is comparable to the strain energy carried by a particle. More specifically, the operator energy in nonlocal operator method describes the irregularity of a function around the crack tip. The irregularity is the part that cannot be described by the continuous function. For continuous domain, the strain energy density is much larger than the operator energy density. However, for particles around the crack tip, the operator energy density is far from zero and the irregularity due to the singularity around the crack tip increases comparably to the strain energy density. In this sense, the operator energy density can be viewed as an indicator for the crack tip.
Unlike the strain energy density, the hourglass energy density describes the irregular deformation around the crack tip. It depends on the penalty for the strain energy. Larger penalty improves the continuity of deformation, but the extent of hourglass energy compared with the strain energy density is hard to estimate. In this paper, we propose a special manner to estimate the critical hourglass strain. Let the critical bond strain be denoted by s max , which may depend on the characteristic length scale of the support, critical energy release rate and the elastic modulus. When the maximal strain reached s max , the damage process is activated and the critical hourglass strain s hg max is set as the maximal hourglass strain s hg ij for all bonds in the computational model. In the sequential calculation, when the hourglass strain of a bond is larger than s hg max , the damage on that bond occurs, which is mathematically described as where d ij denotes the damage status between particle i and particle j.
The damage of a particle is calculated as Every time one particle is removed from the neighbour list, the nonlocal gradient for the central particle should be recalculated based on the remaining "healthy" neighbour. We will apply this rule to model fractures in 2D and 3D linear elastic material.

Numerical implementation
We have applied NOM to derive the nonlocal strong forms for the traditional continuum model in §3. Two representive nonlocal theories, the dual-horizon peridynamics by Eq.44 for fracture modeling and the nonlocal thin plate by Eq.60, are selected for numerical test. For the DH-PD, the focus is on the test of instability criterion for quasi-static fracture modeling by explicit time integration method. The nonlocal thin plate is compared with finite element method. The primary step in the implementation is the calculation of internal force based on the governing equations. In the first step, the computational domain is discretized into particles.
where N is the number of particles in the domain. Then the support of each particle is represented by a list of particle indices, where j is the global index of the particle and n i is the number of particles in S i .
The gradient g ij and Hessian h ij for two particles i, j can be assembled by collecting terms in K i · p ij according to Eq.21 or Eq.23, where with weight function ω(r ij ) = 1/|r ij | 2 . The nonlocal differential derivatives at point i can be calculated as The nonlocal operators in∂u i can be used to define the strain tensor, stress tensor, moment and others.
In discrete form, Eq.44 and Eq.58 become In Eq.93 and Eq.94, the volume of particle i is multiplied on both sides of the equations. It is not required to calculate the internal forces from the dualsupport. Let f i = 0, 1 ≤ i ≤ N denote the initial internal force on particle i. For each particle, one only needs to focus on the support, calculating the forces and adding the force to the particle internal force where a → b denotes the addition of a to b. The process of adding force −ω(r ij 1 )P i ·g ij ∆V j ∆V i to f j is equivalent to accumulating the internal forces from particle j's dual-support.
For the calculating of internal force of thin plate, the same applies In order to maintain the stability of the nonlocal operator, the discrete form of Eq.39 is For particle i with support S i , the hourglass force is calculated as follows When the internal force is attained and the contribution of the external force boundary condition or body force is accumulated, the basic Verlet algorithm [72] outlined as follows is used to update the displacement where u i denotes the displacement or deflection, v i the velocity and a i = f i

Accuracy of nonlocal Hessian operator
We test two cases with analytical function The exact second-order partial derivatives are w ,xx = w ,yy = 2, w ,xy = 0 in 2D (103) w ,xx = w ,yy = w ,zz = 2, w ,xy = w ,yz = w ,xz = 0 in 3D (104)   For regular particle distribution in 2D, small number of particles in support can accurately define the nonlocal Hessian operator, as shown in Fig.2,  Fig.3 and Fig.4. For irregular particle distribution as shown in Fig.5, larger number of particles in support are required to define the nonlocal Hessian operator, as shown in Fig.6 and Fig.7.   For regular particle distribution in 3D, small number of particles in support can accurately define the nonlocal Hessian operator, as shown in Fig.8, Fig.9 and Fig.10.

Single-edge notched tension test
In this example, we tested the nonlocal elasticity by Eq.43 for single-edge notched tension in 2D under plane stress condition. The geometry setup is given in Fig.15. The bottom is fixed while the top of the plate is applied with velocity boundary condition v = 1 m/s, which can achieve the quasistatic condition. The material parameters are E = 210 GPa, ν = 0.3 and critical strain is set as s max = 0.02. The plate is discretized into 100×100 particles. Each particle's support consists of 33 nearest neighbours. The initial crack is created by modifying the neighbour list when searching the    nearest neighbours. The fixed number of neighbours in support results in particles near the boundary with relatively large support sizes and particles in the centre of the plate with small support sizes. A duration of T = 6.5×10 −6 seconds is integrated by approximately 4500 steps at a time increment of ∆t = 1.5418 × 10 −9 seconds. The displacement field u y at step 3250 and step 4200 are depicted in Fig.16(a) and Fig.16(b), respectively.   Fig.17(a) is the displacement field u y at full damage, where the interaction of internal force between the two half planes is cut and rigid body displacement dominates. Fig.17(b) is the distribution of hourglass energy. We can observe that the hourglass energy is concentrated on the crack surface and crack front tip. Fig.17(c) and Fig.17(d) are the snapshots of damage field, which confirms that the instability criterion in §4 is stable for fracture modelling. Although the plate is solved by an explicit dynamic method, the kinetic energy is much lower than the strain energy as shown by Fig.18(a). The dynamic load curve agrees well with that by the finite element method in Ref [70], as shown by Fig.18(b). One possible reason for the difference in reaction force increment is due to explicit algorithm and nonlocal effect of current formulation.

Out-of-plane shear fracture in 3D
For brittle fracture, the basic modes of fracture are tensile fracture, in-plane shear fracture, out-of-plane shear fracture. In this section, we apply the instability damage criterion to the out-of-plane shear fracture, as shown in Fig.19. The dimensions of the specimen are 5×2×1 mm 3 , as shown in Fig.20. The size of the initial crack surface is 2 × 1 mm 2 . The velocity boundary conditions u z = 1 m/s are applied. The model is discretized into 86961 particles with particle size ∆x = 0.05 mm. Each particle has 102 neighbours in its support. Material parameters include elastic modulus E = 210 × 10 9 Pa and Poisson ratio ν = 0.3 and density ρ = 7800 kg/m 3 . The time step is selected as ∆t = 7.7 × 10 −9 seconds. A total of 3000 steps are calculated. The crack surface starts to propagate at step 1550. The crack surface at different steps are depicted in Fig.21.

Conclusion
In this paper, we employ the recent proposed NOM to derive the nonlocal strong forms for various physical models, including elasticity, thin plate, gradient elasticity, electro-magneto-elastic coupled model and phase field fracture model. These models require second order partial derivative at most    and we make use of the second-order NOM scheme, which contains the nonlocal gradient and nonlocal Hessian operator. Considering the fact that most physical models are compatible with the variational principle/weighted residual method, we start from the energy form/weak form of the problem, by inserting the nonlocal expression of the gradient/Hessian operator into the weak form, based on the dual property of the dual-support in NOM, the nonlocal strong form is obtained with ease. Such a process can be extended to many other physical problems in other fields. The derived strong forms are variationally consistent and allow elegant description for inhomogeneous nonlocality in both theoretical derivation and numerical implementation. We also propose an instability criterion in nonlocal elasticity or dualhorizon state-based peridynamics for the fracture modeling. The criterion is formulated as the functional of nonlocal gradient in support, which minimizes the zero-energy deformation that cannot be described by the nonlocal gradient. Such an operator functional approaches zero for continuous fields but has comparable value to the strain energy density for the deformation around the crack tip. During the fracture modeling by removing particles from the neighbor list, it is safer to delete the particle with larger zero-energy deformation. The numerical examples for 2D/3D fracture modeling confirm the feasibility and robustness of this criterion. The instability criterion is possible applicable for anisotropic elastic material and hyperelastic materials.
A A simple example to illustrate dual-support In order to facilitate the comprehension of dual-support, let us consider 4 particles in Fig.22, each with particle volume ∆V i , i = {1, 2, 3, 4} and Ω = 4 i=1 ∆V i . Obviously, the support and dual-support can be listed as follows. Here we neglect whether the shape tensor is invertible or not. The most common formula in the derivation based on NOM and variational principle is the double integrations in support and whole domain.

Consider the double integrations
At last, we obtain Above equation is widely used in the derivation of nonlocal strong form from weak form. Such expression is valid in the continuum form as well as in discrete form.