Mixed peridynamic formulations for compressible and incompressible ﬁnite deformations

The large ﬂexibility of meshfree solution schemes makes them attractive for many kinds of engineering applications, like Additive Manufacturing or cutting processes. While numerous meshfree methods were developed over the years, the accuracy and robustness are still challenging and critical issues. Stabilization techniques of various kinds are typically used to overcome these problems, but often require the tuning of unphysical parameters. The Peridynamic Petrov–Galerkin method is a generalization of the peridynamic theory of correspondence materials and offers a stable and robust alternative. In this work, the stabilization free approach is extended to three dimensional problems of ﬁnite elasticity. Locking-free mixed formulations for nearly incompressible and incompressible materials are developed and investigated in convergence studies. In general, an efﬁcient implicit quasi-static framework based on Automatic Differentiation is presented. The numerical examples highlight the convergence properties and robustness of the proposed formulations.


Introduction
In last decades many different methods for the approximate solution of partial differential equations were developed and applied to various engineering problems. The Finite Element Method (FEM) is well accepted for the solution of a wide range of problems in industry. However, for complex evolving domains it is advantageous to have a more flexible discretization scheme. For instance in the fast growing field of Additive Manufacturing like Selective Laser Melting, the production process is accompanied by evolving surfaces, large deformations and phase changes. A promising discretization approach consists in meshfree methods which include the Smoothed Particle Hydrodynamics (SPH) introduced by Lucy [24] and Gingold and Monaghan [13], the Reproducing Kernel Particle Method (Liu et al. [23]), the Element Free Galerkin method (Belytschko et al. [2]) and the more recent non-ordinary state-based peridynamic correspondence approach (see Silling et al. [31]) where B T. Bode bode@ikm.uni-hannover.de 1 Institute of Continuum Mechanics, Leibniz University Hannover, Appelstr. 11, 30167 Hannover, Germany each of them has their specific advantages and drawbacks. Peridynamics is a non-local field theory based on integrodifferential equations which is idealized and widely used to model discontinuities as occurring in fracture simulations (see e.g. Silling [28] and Silling [29]) and in recent years the area of applications rapidly increased. For a review on Peridynamics the reader is referred to Javili et al. [17].
Besides the advantage of flexibility in discretization, the efficiency and accuracy of the meshbased FEM is hard to achieve. Some reasons are consistency, stability and the imposition of boundary conditions which are challenging aspects in the scope of meshfree methods. Rank deficiency which result in spurious zero or low energy modes is often addressed by corrections to gain stability (see e.g. the works of Littlewood [22], Breitenfeld et al. [5], Silling [30], Weißenfels and Wriggers [35]). However, unphysical parameters have to be determined that can be case sensitive and reduce the accuracy of the solution.
In the context of the peridynamic correspondence model, approaches were developed to overcome the rank deficiency without the use of such corrections in the last years.
Tupek and Radovitzky [33] proposed an extended correspondence formulation based on Seth-Hill strains. A higher order approximation using modified weight functions has been developed in Yaghoobi and Chorzepa [38]. Chowdhury et al. [11] adressed the low-energy modes by means of separation of the so-called family which is the neighborhood of a particle (see Fig. 1) into subdivisions. The use of bond-level and bond-associated deformation gradients is studied in Breitzman and Dayal [6], Chen [9], Gu et al. [14], Madenci et al. [27] and Chen and Spencer [10]. A stresspoint method to overcome the rank deficiency was presented in Luo and Sundararaghavan [25]. Hillman et al. [15] unified the peridynamic deformation gradient with the implicit gradient approximation, studied the convergence of the differential operations and proposed a reproducing kernel peridynamic method.
In a further approach the Peridynamic Petrov-Galerkin (PPG) method was introduced, see Bode et al. [4]. It is based on the principle of virtual displacements applied to the peridynamic momentum equation and depicts good accuracy and suffers no oscillations when applied within large deformations. In case of linear ansatz functions it reduces to the commonly used correspondence formulation of Peridynamics. Hence it is applicable to any local material models. In addition, the Petrov-Galerkin formulation enables to a certain extent knowledge transfer from the well established FEM.
In this paper the PPG method is extended to three dimensions and more efficient interpolating Moving Least Square shape functions are utilized. Furthermore, the regime of incompressible material behavior is investigated. Like in the FEM, pure displacement-based peridynamic approaches exhibit locking phenomenas as e.g. a stiffer response or divergence. Mixed displacement-pressure formulations are known to overcome these problems and possess good convergence rates (see Simo et al. [32]). Therefore, mixed PPG formulations for weakly compressible and incompressible material behavior are presented and the convergence and robustness is compared to Finite Elements. The resulting formulations can be viewed as a non-ordinary state-based correspondence model incorporating local constraint equations by means of Lagrange multipliers. The structure of the paper is as follows. The governing equations for nonlinear elasticity are presented in Sect. 2. Section 3 presents the PPG discretization and formulates the mixed formulations based on Automatic Differentiation (AD). Numerical convergence studies for compressible and incompressible benchmark test as well as a numerical inf-sup test are considered in Sect. 4. The paper closes with concluding remarks in Sect. 5.

Governing equations for finite elasticity
The peridynamic correspondence theory is based on the nonlocal peridynamic momentum equation. However, it can be applied to constitutive models from the local theory. There-fore, the local kinematic description of the continuum can be used (see e.g. Wriggers [36] and Holzapfel [16]). For an elastic body the current position of a material point x is given by with the initial position X and the displacement u. The deformation gradient is defined as where the Jacobian J describes the volumetric part. We also define the right Cauchy-Green tensor C by and its isochoric part as The body is subject to the momentum equation at each particle X. It states in the non-local peridynamic form as an integral over momentum transfers with neighboring material points. Herein, t and t are called the pairwise force densities. While t stands for the force acting on X exerted by a neighboring material point X inside the family H (see Fig. 1), the pairwise part t arises from the collective deformation of family H of X and is considered by means of Newton's third law. ρ 0 b is a body force acting on X. The local counterpart is defined with respect to the initial configuration as Div P + ρ 0 b = 0 (6) with the first Piola-Kirchhoff stress tensor P. A correspondence formulation t = t (P) links the first Piola-Kirchhoff stress with a state of pairwise force densities in the discretized form (see for instance Silling et al. [31], Madenci and Oterkus [26] and Bode et al. [4]). By applying the principle of virtual displacements on equations (5) and (6) the peridynamic virtual strain energy yields and the local counterpart excluding surface forces Using the correspondence formulation, the constitutive laws can be based on the local theory for elastic finite deformations. In case of incompressible material behavior the Neo-Hookean strain energy function is used and in case of compressible material behavior extended to The bulk modulus K and the shear modulus μ can be calculated from the Young's modulus E and the Poisson's ratio ν or lame constants λ and μ. The derivation of the strain energy function with respect to the deformation gradient yields the first Piola-Kirchhoff stress

Formulation of the peridynamic Petrov-Galerkin method
The PPG method is based on a subdivision of the body into a finite number of discrete particles. Considering collocation, these particles serve both as material points and nodes (see Fig. 1). The connectivity between particles is established by so-called families which define the surrounding neighborhood of each particle. A particle k interacts only with the neighbors inside its family H k . The local linear momentum preservation is ensured if these interactions are always pairwise, i.e. each particle j inside family H k has to have particle k in its own family H j . Such a pairwisity is usually fulfilled by having spherical families with a radius of the so-called horizon δ which is constant inside body . However, the PPG method is not restrictive to spherical family shapes. In this work, for irregular particle distributions the families are made up of the nearest N = 20 particles (for three dimensions). To ensure pairwisity, all pairwise particles which are not already inside a family are added afterwards. In case of uniform particle patterns, a horizon δ of 1.51 times the particle spacing is used. The non-local kinematics of a family H k consists of the bonds X k j = X j −X k , i.e. the distance vector from particle k to its neighbor j. Analogously, the deformed bonds x k j = x j − x k and the displacement differences u k j = u j − u k = x k j − X k j are defined. For enforcing local constitutive laws, the non-local kinematics have to be transformed to a local defined measure like the deformation gradient F. This transformation is performed by the following approximations for the actual and virtual deformation gradient: Fig. 1 The connectivity between particles with their representative volume V is established by the family H and the bonds X. Internal forces acting on a particle are considered in terms of the pairwise force densities t where F k j stands for the deformation gradient at the position of particle j with respect to family H k . As the virtual deformation gradient is evaluated only at the master particle k a single superscript is used. The i-th ansatz function of family H k at the position of particle j is written as N k ji , where the subscript u and δu stand for the trial and test shape function, respectively. With this ansatz, the aforementioned correspondence formulation can be derived. It relates the first Piola-Kirchhoff stress tensor P k j from a local material model to the non-local pairwise force densities t k j . This relation is based on the equivalence of the virtual work of Eq. (7) and (8) and states Inserting the pairwise force density relation into the peridynamic momentum equation and integrating over the particle volume V k yields the discretized form in the PPG method at particle level

PPG shape functions
The choice of shape functions for test and trial functions is a crucial point in the PPG method. It can decide over rank deficiency, stability, the accurate imposition of Dirichlet and Neumann boundaries and the convergence rate. In Bode et al. [4] it was shown that in case of a linear Weight Least Square (WLS) functions, the PPG method reduces to the nonordinary state based peridynamic method for correspondence materials.
The key features for convergence to an accurate solution consists of consistency, i.e. the property to exactly reproduce a polynomial space (see Krongauz and Belytschko [20] and Belytschko et al. [3]). Further necessary requirements on discretization schemes consist in the absence of rank deficiency and an exact imposition of boundary conditions. While at least linear consistency is fulfilled by most meshfree shape functions, the second and third aspects are more challenging for meshfree methods. Rank deficiency and the resulting low-energy modes are oftentimes addressed by stabilization techniques which usually include a tuning parameter and lead to an artificial stiffness. In general the PPG method does not need a stabilization, provided that the trial function ansatz is not linear which would lead to a constant deformation gradient inside a family. This is the case for many shape functions as e.g. second order WLS functions, Moving Least Square (MLS) functions, Local Maximum Entropy (LME) approximants (see Arroyo and Ortiz [1]).
Regarding the exact imposition of boundary conditions, one has to differentiate between Dirichlet and Neumann boundaries. The latter one can be imposed accurately by fulfilling the integration constraint (see Bode et al. [4] and Weißenfels [34]). The peridynamic form states where A k is the surface of particle k. This value is only non zero, if the particle is located at the surface of the body . The vector N k is the corresponding unit normal vector at the boundary node. Dirichlet boundaries can be imposed accurately by using the Kronecker-δ property of the test functions (see also Bode et al. [4]). In the following, an interpolating MLS approach based on differences is proposed.
WLS and MLS are fitting techniques for arbitrary point clouds (see Lancaster and Salkauskas [21]). Within these techniques, a field is approximated by a linear combination of shape functions, like for instance monomials. The coefficients are determined in such a way that the squared error at the data points is minimized. In case of WLS the errors are weighted statically and in case of MLS dependent on the evaluation point. For the evaluation at particle j inside family H k this is performed via the minimization of the functional where the approximation errors ki = û k X k j − u ki , volume weights V i and distance weights ω ji were introduced. As in the PPG method shape functions are only used to approximate derivatives, the fitting can be done with respect to differences which leads to an enforcement of the interpolation at the reference particle. In a MLS sense, this perspective X j can move with the according weights ω ji which results in the ansatz monomials of first order The distance weights ω ji regulate the influence of the neighboring particles. A localized fitting can be introduced by a weighting function that decreases with increasing distance. The weights can be set to with p = 2, or in case of regular meshes even more localized with p = 5. From Eq. (16), applied to the displacement difference field in family H k , follows the approximation where N k ji is the i-th ansatz function of family H k evaluated at particle j. The mass matrix, which is also called shape tensor in Peridynamics or correction matrix in SPH, ensures consistency up to the order of completeness of the ansatz monomials. It depends on the position inside the family and is given by Utilizing zeroth order consistency, to the gradient of the displacement difference approximation (19) yields the standard approach for the displacement gradient with respect to the reference configuration The gradient of the shape functions ∂ N k ji ∂X at data points follows from (19) and the application of the zeroth order consistency condition

Mixed PPG formulations
By integrating the peridynamic virtual strain energy (7) over the body , a general weak form can be formulated: With this an implicit quasi-static framework can be derived from the discretized form by separating the body into a set of n p particles with the corresponding volume of V k . Then, the approximate weak form yields The discretized virtual strain energy δW k is based on the concept of pairwisity exploiting the correspondence formulation of Eq. (13) Utilizing the extended strain energy function (10) yields a displacement based PPG formulation. In case of using MLS trial and test functions of first order, this formulation will be termed "PPG Unl", since the approach for the actual displacements u is nonlinear (nl).
Next, a mixed displacement-pressure-dilation formulation motivated by the Q1P0/H1P0 finite elements (see Simo et al. [32]) is presented. It is based on an additive split of the strain energy into an isochoric and volumetric part which leads to an additive split of the pairwise force densities The isochoric part is handled analogously to the pure displacement based formulation using the isochoric Neo-Hookean strain energy function of Eq. (9). The volumetric part is now derived from the Hu-Washizu potential By applying particle discretization the discrete form states where the pressure p k and the dilation k are independent variables of particle k. The resulting formulation will further be denoted as "PPG UnlP0". In case of incompressible material behavior, i.e. ν = 0.5, the dilation k = 1 which yields the "PPG IUnlP0" formulation. Both formulations satisfy the conservation of linear momentum as the volumetric residual results in pairwise force densities t k j vol . A proof is provided in "Appendix A", which also shows the additive split of the pairwise force densities into

Residual and tangent
The weak form (24) is minimized by means of the Newton-Raphson method. Therefore, the global residual vector and the tangent matrix can be either computed numerically (see e.g. Brothers et al. [7]) or by a consistent linearization (see Bode et al. [4]). In this work, the residual and tangent matrix are assembled family-wise where the linearization is performed with the symbolic Automatic Differentiation (AD) tool AceGen, see Korelc and Wriggers [19]. Hence, for the PPG Unl formulation the residual and tangent matrix follows Note, that the first derivation of the weak form leads back to the residual of Eq. (14). In addition to the family-wise residual and stiffness matrix (30), the part due to volumetric forces for the mixed PPG UnlP0 and PPG IUnlP0 formulations can again be derived with AD to UnlP0: Mixed displacementdilation-pressure approach; Constant pressure and dilation Weight function ω 2 (regular and distorted), ω 5 (regular) (18) Therefore, the family-wise stiffness matrix for the PPG UnlP0 and PPG IUnlP0 approach yields In the first case on the left of Eq. (32), the pressure and the dilation can be eliminated at family-level by means of static condensation which leads to a higher efficiency and a better conditioning.

Numerical examples
In this section, the new formulations are tested by means of several examples. First, the locking behavior is investigated. Second, the convergence is compared to FEM solutions for regular and irregular particle distributions. Third, the robustness is demonstrated in a torsion test and finally, a numerical inf-sup test is performed.

Punch problem
The first example investigates a punch into a block of solid material as shown in Fig. 2 (cf. Wriggers [37]). The block has the dimensions 0. Free moving particles are colored green. When applying the load, the block undergoes large deformations. In Fig. 2 the final configuration as well as a cross section is depicted for the PPG IUnlP0 formulation with a discretization of 24948 particles. The displacement in z-direction is colored, see legend in Fig. 2d. The convergence of the nodal displacement in z-direction in the middle of the block at the upper surface is analyzed. Figure 3 depicts the convergence behavior for the different formulations using MLS1 with ω 5 as both test and trial functions, see Table 1. The displacement based PPG Unl formulation is only pictured for a Poisson's ratio of ν = 0.4 as it fails to converge or does almost not deform due to locking in the incompressible regime. However, the PPG UnlP0 formulation as well as the incompressible PPG IUnlP0 formulation are free of locking and converge at about the same rate as those of FEM H1P0.

Cook's membrane problem
The second example is Cook's membrane problem. As depicted in Fig. 4, a tapered cantilever beam is clamped on the left (fixed purple particles) and loaded at the right (red particles) in a single step with a load of q = 4 N m 2 in positive z-direction. The block is modeled as a Neo-Hookean solid with the Lame constants μ = 40 N m 2 and λ = 100 N m 2 . The final configuration in case of the PPG UnlP0 formulation with a regular discretization using 95, 904 particles is shown in Fig. 4 where the displacement in z-direction is colored.
In a convergence study the PPG UnlP0 formulation compares regular meshes, using ω 2 and ω 5 distance weights, with Different discretizations and weight functions in the PPG method are compared to the mixed H1P0 Finite Elements an irregular particle distribution (PPG UnlP0-d) using the ω 2 weight function. The finest irregular and regular particle distribution are depicted in Fig. 5 showing the von Mises stress in the final configuration. The stresses are highly concentrated at the weak singularity at the upper left edge and exhibit a smooth transition over the whole beam. Figure 6 shows the displacement of the right upper front edge over the number of nodes. The convergence rate of the PPG formulations is comparable to the H1P0 Finite Element while it converges from the opposite side. In case of strongly localized ω 5 distance the converged displacement matches the one of FEM, for more smoothing ω 2 weights a slight discrepance is observed.

Torsion problem
The third example demonstrates the robustness of the mixed PPG UnlP0 in a torsion test given in Kadapa et al. [18]. The Via adaptive load stepping, the twisting angle is increased until the convergence fails for a further increment of at least one degree. The total angle as well as the maximal angle that can be applied in a single load step are shown in Table 2 for a discretization into 16 × 16 × 5 · 16 particles and elements, respectively. Both measures show the superiority of the mixed formulations over the pure displacement based approaches. While the FEM H1 and PPG Unl formulations need small load steps which soon fall below one degree, the FEM H1P0 and especially the PPG UnlP0 approaches can deal with very large load steps. The peridynamic approach is even more robust and effective for large deformations compared to standard finite element methods.

Numerical inf-sup test
In mixed FE methods, the inf-sup or Ladyzhenskaya-Babuška-Brezzi condition ensures together with the ellipticity condition the existence, uniqueness and stability of the regarding discretization. For nearly incompressible material behavior in linear elasticity, the inf-sup condition (see Fortin and Brezzi [12]) states  where P h and U h are finite dimensional spaces of trial pressure and displacement fields, respectively. Due to the complexity of meshless shape functions, an analytical proof whether the condition is satisfied is difficult to state. However, Chapelle and Bathe [8] proposed a numerical test which can give a prediction for the satisfaction. Therefore, a bloc of material with the applied essential boundaries as shown in Fig. 9a is discretized by a series of regular and distorted particle distributions consisting of N × N particles or elements with N = {2, 4, 8, 16, 32} (see Fig. 9b, c). The limit β of the inf-sup condition of Eq. (33) can be determined from the following generalized eigenvalue problem: where U h is the global vector of nodal displacements. The matrix S h results from the H 1 norm of the discrete displacement field and T h corresponds to the discrete pressure projection Considering the range of small strains and using static condensation, it states The inf-sup value of a particular resolution computes to the smallest non-zero eigenvalue of the generalized eigenvalue problem defined in Eq. (34): where k − 1 is the number of zero-eigenvalues. The number of spurious pressure modes can be determined from with the number of global displacement and pressure degrees of freedom n u and n p . Figure 8 depicts the inf-sup values of  Fig. 9, where the logarithm of the inf-sup value is plotted over the logarithm of the inverse of the number of particles in each dimension for regular and distorted (-d) particle and node distributions, respectively. While the Q1P0 element exhibits an inf-sup value converging towards zero, the inf-sup value of the PPG UnlP0 formulation is bounded from below by a non-zero value the model problem of Fig. 9. Whereas the Q1P0 element violates the inf-sup condition, the UnlP0 formulation passes the numerical inf-sup test as the inf-sup value is bounded from below and k pm = 0 for both uniform and irregular particle distributions. In case of the distorted particle distribution, the same neighborhood as in the according regular distribution is choosen. The satisfaction of the numerical inf-sup test is in agreement with the presented examples where no locking phenomenas in the range of nearly incompressible material behavior were observed.

Summary
In this work, the novel Peridynamic Petrov-Galerkin method is extended to three dimensional problems and a special kind of interpolating Moving Least Square approximation of derivatives based on differences is proposed. The locking behavior of a pure displacement based formulation is investigated and by transferring knowledge from the Finite Element technology mixed displacement-dilation-pressure formulations are developed. An implicit quasi-static framework based on efficient Automatic Differentiation is presented.
The numerical examples show that the proposed approach achieves good performance concerning robustness and convergence while maintaining the advantages of meshfree particle methods.
Acknowledgements Open Access funding provided by Projekt DEAL.
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://creativecomm ons.org/licenses/by/4.0/.

A Proof of momentum preservation
The preservation of linear momentum is ensured by the concept of pairwise forces. As the isochoric part of the strain energy is handled in the standard pure displacement based way, the momentum preservation is automatically fulfilled and it remains to show that the volumetric part of the residual can also be written in terms of pairwise force densities. Thus, starting from the local volumetric potential On the left: Model problem for numerical inf-sup test. On the right: exemplary discretizations for regular and irregular particle distributions. The Dirichlet boundary conditions are applied by and additional layer of wall particles. The displacements of blue colored particles are fixed horizontally and the purple particle is fixed in both horizontal and vertical direction the residual contribution to particle j can be computed as the derivative with respect to the particles displacement ∂U k vol ∂u j = V k p k ∂ J k ∂F k : where the chain rule is used. Inserting the deformation gradient approach (12) and switching to index notation yields The resulting pull-back of the volumetric Cauchy-stress can be exchanged by the first Piola-Kirchhoff stress tensor.
Analogous to the contribution to the nodal residual of particle j, the corresponding contribution to particle k follows where the zeroth order consistency conditions is utilized. The internal force acting on particle k can now be assembled from the potentials of its neighboring particles and rearranged to where the volumetric pairwise force density states As the inner forces acting on a particle can be expressed with pairwise force densities, the linear momentum preservation holds if pairwisity is fulfilled.