Derivation of dual-horizon state-based peridynamics formulation based on Euler–Lagrange equation

The numerical solution of peridynamics equations is usually done by using uniform spatial discretisation. Although implementation of uniform discretisation is straightforward, it can increase computational time significantly for certain problems. Instead, non-uniform discretisation can be utilised and different discretisation sizes can be used at different parts of the solution domain. Moreover, the peridynamic length scale parameter, horizon, can also vary throughout the solution domain. Such a scenario requires extra attention since conservation laws must be satisfied. To deal with these issues, dual-horizon peridynamics was introduced so that both non-uniform discretisation and variable horizon sizes can be utilised. In this study, dual-horizon peridynamics formulation is derived by using Euler–Lagrange equation for state-based peridynamics. Moreover, application of boundary conditions and determination of surface correction factors are also explained. Finally, the current formulation is verified by considering two benchmark problems including plate under tension and vibration of a plate.

Analytical solutions of PD equations are not usually possible. Therefore, numerical techniques, especially meshless method, are widely utilised. In majority of the publications in the literature, uniform discretisation is used for spatial discretisation. Although uniform discretisation is simple to implement, for certain applications, it unnecessarily increases computational time since only some parts of the solution domain require a finer discretisation and other parts can be discretised by using a coarser grid. Hence, non-uniform discretisation can be very beneficial in terms of computational time. In addition to non-uniform discretisation, horizon size can also be different at different parts of the solution domain either to reduce the computational time or to capture correct physics of the problem. However, this requires extra care especially to satisfy conservation laws. To accomplish these requirements, Ren et. al. [45,46] introduced dual-horizon peridynamics concept. In this study, we present a different way to derive dual-horizon peridynamics formulation for state-based peridynamics by utilising Euler-Lagrange equations. Moreover, application of boundary conditions and determination of surface correction factors are also explained. Finally, the current formulation is verified by considering two benchmark problems including plate under tension and vibration of a plate.

Dual-horizon peridynamics formulation based on Euler-Lagrange equation
In this study, dual-horizon peridynamics formulation is derived based on Euler-Lagrange equation. Peridynamics is a non-local continuum mechanics formulation. Therefore, the body is composed of small volumes called "material points" as shown in Fig. 1.
The equation of motion of a particular material point k at location x (k) in the initial configuration can be obtained by using Euler-Lagrange equation as [4] d dt where L is the Lagrangian which represents the difference between the total kinetic energy of the body and the total potential energy of the body, t represents time, u (k) andu (k) are the displacement and velocity of the material point k. The total kinetic energy of the body can be calculated by summing the kinetic energies of all Fig. 1 Material points and the horizon concept [4] material points in the body as where ρ is the density and V (i) is the volume of the material point i. On the other hand, the total potential energy of the body is the difference between the total strain energy of the body and the energy due to external forces which can be written as where W is the strain energy density which can be defined for the material point k as with ω being the micropotential arising due to interaction between material points k and j, N k and N j represent the number of family members within the horizon of material points k and j, respectively, and y is the position of the material point in the deformed configuration. Note that micropotentials ω (k)( j) and ω ( j)(k) do not need to be equal to each other since they have a different domain of influence, horizon (Fig. 1). In Eq. (3), b is the body load acting on a particular material point. Substituting Eqs. (2) and (3) in Euler-Lagrange equation given in Eq. (1), the equation of motion of the material point kcan be obtained as where peridynamic forces given in Eq. (5) can be written in terms of micropotentials between interacting material points, ω, as For ordinary state-based peridynamics, peridynamic forces can be rewritten for variable size horizon case as with α k j = 1, ω k j = 0 0, ω k j = 0 (10) The peridynamic dilatation in Eqs. (8,9) for the material points k and j is defined as The parameters a, b and d are peridynamic parameters which can be expressed in terms of material constants of CCM for an isotropic material as with κ and μ being the bulk and shear moduli, respectively, h is the thickness of the two-dimensional domain and δ is the horizon size of a particular material point. Note that in Eqs. (10) and (11), the micropotential being zero means one of the material points is not within the horizon of the other material point. As shown in Fig. 2, although the material point with a blue (smaller) horizon is inside the horizon of the material point with red (larger) horizon, the material point with a red horizon is outside of the horizon of the material point with a blue horizon.
For non-ordinary state-based peridynamics, peridynamic forces for variable horizon can be expressed as where P is the first Piola-Kirchhoff stress tensor. The shape tensor, K, for the material points k and j is defined as where the symbol ⊗ represents the dyadic product. Note that non-ordinary state-based formulation can suffer from the zero-energy mode problem. In this study, the zero-energy mode problem is prevented by following the approach given in Silling [47]. For bond-based peridynamics, peridynamic force between two interacting material points is only influenced by the motion of these two material points. In other words, only micropotentials which belong to both interacting points exist, i.e. ω ik = ω ki = 0, if i = j. Therefore, peridynamic forces given in Eqs. (6) and (7) can be simplified as These forces can also be expressed as where c is the bond constant which can be written for an isotropic material and two-dimensional domain as with E being the elastic modulus.

Boundary condition
Since PD equations are in the form of integro-differential equations, the way to apply boundary conditions in PD is different than CCM. Boundary conditions in PD are applied through volumes, and these volumes are sometimes created as fictitious regions outside the actual solution domain. In this study, boundary volumes are chosen as fictitious regions, f , outside of the main solution domain, , as shown in Fig. 3, which is an effective procedure suggested in [14].

Displacement constraints
In this study, the size of the fictitious region is specified as twice the horizon size, i.e. 2δ, and only twodimensional problems are considered. However, the presented approach can be easily extended to threedimensional models. The prescribed boundary value of the displacements U * and V * in the x− and y− directions can be applied by specifying the displacements of the material points in the fictitious region, u f and v f , in terms of displacements of the material points in the actual domain, u and v, as (Fig. 3)

Traction boundary conditions
Traction boundary conditions can also be applied similar to displacement constraints by using a fictitious region as shown in Fig. 4. The displacements of the material points in the fictitious region depend on the unit normal of the boundary on which the traction boundary condition, σ * i j , with i, j = x, y, is exerted. For the traction boundary with a unit normal in the x−direction, they can be written as [14] u Similarly, for the traction boundary with a unit normal in the y−direction, they can be written as

Surface correction
Determination of surface correction factors for ordinary state-based peridynamics is important not only for material points close to the surfaces due to lack of interactions, but also for material points inside the solution domain since spatial numerical integration scheme based on meshless approach introduces numerical error to the solution. These errors can be corrected by using surface correction factors and modifying the peridynamic dilatation terms and peridynamic force expressions given in Eqs. (8,9,13,14) as and where β and λ represent the surface correction terms for the dilatation and peridynamic force expressions, respectively. As explained in Madenci and Oterkus [4], surface correction factors can be obtained by considering a simple loading condition and determining the ratio of the dilatation and strain energy density of a particular material point calculated from CCM and PD. Note that surface correction factors are not necessary for the non-ordinary state-based formulation since correction is automatically done within the formulation.

Numerical results
As mentioned earlier, although uniform discretisation is widely used in the peridynamic studies available in the scientific literature, it can be computationally advantageous to have flexibility to use different grid sizes at different parts of the solution domain which is a common procedure used in other numerical methods such as finite element method (FEM). Hence, in this section, we present the capability of dual-horizon peridynamics

Plate under tension
In the first problem case, a square plate with dimensions of L = W = 1m is considered. The plate has a thickness of 0.01 m and subjected to a uniaxial tension loading of σ * = 200 MPa at the right edge as shown in Fig. 5. The applied loading is specified by creating a fictitious region at the right edge as demonstrated in Fig. 6. Elastic modulus and density of the plate are specified as 200 GPa and 7850 kg/m 3 , respectively. The left edge of the plate is fully fixed by using a fictitious region (Fig. 6). The steady-state solution is obtained by utilising the adaptive dynamic relaxation technique [48].

Ordinary state-based peridynamics
A plate under tension is first analysed by using ordinary state-based peridynamics. The solution domain is split into two equal regions as Region 1 and Region 2 as shown in Fig. 7. The discretisation size in Region 1, In addition, the effect of the horizon size is investigated by considering the same horizon sizes for Regions 1 and 2 as shown in Fig. 10, i.e. δ 1 = δ 2 = 6 x 1 = 3 x 2 , and compared against the different horizon case considered earlier with the horizon size in Region 1, δ 1 = 3 x 1 , and the horizon size in Region 2, δ 2 = 3 x 2 . As shown in Figs. 11 and 12, although both cases provide good agreement with FEM results for horizontal displacements, equal horizon case yields better vertical displacement results with respect to FEM results since there are more material points inside the horizon in Region 1 for this case.

Non-ordinary state-based peridynamics
In this section, a plate under tension problem is analysed by using non-ordinary state-based peridynamics. As in the previous section, the horizon size and discretisation size of Region 1 are half of the Region 2 (Fig. 13). Five different horizon size cases are considered, i.e. δ 1 = n x 1 , δ 2 = n x 2 , n =1,…5, and the variation of horizontal and vertical displacements along the central axes is shown in Figs. 14 and 15. According to these figures, the horizon size value of δ 1 = 2 x 1 & δ 2 = 2 x 2 has better agreement with FEM results since it provides smoother transition between the two regions which is a different outcome than ordinary state-based case. As the horizon size increases, error at the interface region increases. Next, equal size horizon for both Regions 1 and 2 is considered and compared against the case with different horizon sizes in Regions 1 and 2 where the horizon size in Region 1, δ 1 = 2 x 1 , is half of the horizon size in Region 2, δ 2 = 2 x 2 (Fig. 16). As shown in Figs. 17 and 18, both cases yield good agreement with FEM results.

Vibration of a plate
In the second problem case, the square plate considered in the previous case is subjected to an initial condition in the form of uniaxial strain of 0.001 in the x-direction as shown in Fig. 19. The solution is obtained using explicit time integration [4] with a time step size of 1 × 10 −7 s (Fig. 20).

Ordinary state-based peridynamics
As in the plate under tension problem, first solution is obtained by using ordinary state-based peridynamics. The solution domain is split into two equal regions as Region 1 and Region 2 as shown in Fig. 21. The discretisation size in Region 1, x 1 = 0.005 m, is half of the discretisation size of Region 2, x 2 = 0.01 m. Moreover, the horizon size of material points in Region 1, δ 1 , is half of the horizon size in Region 2, δ 2 . Next, the effect of discretisation size is investigated. By keeping the same horizon sizes as before, the discretisation size in Region 1 is specified as half of the previous case as shown in Fig. 24. Therefore, the grid size in Region 2 is x 2 = 0.01, whereas the grid size in Region 1 is x 1 = 0.005 for Case 1 and x 1 = 0.025 for Case 2. The horizon size in Region 1 is half of the horizon size in Region 2, δ 2 = 3 x 2 . As shown in Figs. 25 and 26, although horizontal displacements in both cases are similar to each other, Case 2 yields better vertical displacement results with respect to Case 1 since there are more material points inside the horizon of Region 1 in Case 2.

Non-ordinary state-based peridynamics
Finally, the vibration of a plate problem is analysed by using non-ordinary state-based peridynamics. Five different horizon sizes are considered, δ 1 = n x 1 , δ 2 = n x 2 , n =1,…5. The horizon size in Region 1 is half of the horizon size in Region 2 (Fig. 27). As shown in Figs. 28 and 29, the horizon size value of δ 1 = 2 x 1 & δ 2 = 2 x 2 has better agreement with FEM results.
Next, two different non-uniform discretisation cases are compared (Fig. 30). The grid size in Region 2 is x 2 = 0.01, whereas the grid size in Region 2 is x 1 = 0.005 for Case 1 and x 1 = 0.025 for Case 2. The horizon size in Region 1 is half of the horizon size in Region 2, δ 2 = 2 x 2 . As shown in Figs. 31 and 32, both cases provide similar results and good agreement with FEM solution.

Conclusions
In this study, dual-horizon peridynamics formulation was derived by using Euler-Lagrange equation. Dualhorizon peridynamics is an effective methodology when non-uniform discretisation and variable horizon sizes are considered either due to computational efficiency purposes and/or due to physical requirement of the problem. After deriving the dual-horizon peridynamic formulation for both ordinary and non-ordinary state-based peridynamics, two benchmark problems were analysed to investigate the effect of discretisation size and horizon size. By comparing PD results against FEM results, non-ordinary state-based peridynamics provides better agreement for smaller horizons sizes and error at the interface becomes larger as the horizon size increases. On the other hand, for ordinary state-based peridynamics results, the results at the interface are smoother and the results are getting better as the horizon size increases. Moreover, as the number of points inside the horizon is increasing, a better agreement is obtained for a constant horizon size in the solution domain. Finally, with the advancement of additive manufacturing technologies, fabrication of complex materials such as microstructured materials is possible. PD theory can be a suitable alternative for the analysis of microstructured materials to some other approaches presented in the literature [49][50][51], especially for damage prediction.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-18-1-7004.