Peridynamic Petrov–Galerkin method: A generalization of the peridynamic theory of correspondence materials

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 358 (2020) 112636 www.elsevier.com/locate/cma

Peridynamic Petrov–Galerkin method: A generalization of the peridynamic theory of correspondence materials T. Bode ∗, C. Weißenfels, P. Wriggers Institute of Continuum Mechanics, Leibniz University Hannover, Appelstr. 11, 30167 Hannover, Germany Received 14 May 2019; received in revised form 11 September 2019; accepted 12 September 2019 Available online xxxx

Highlights • An ansatz function based correspondence formulation eliminates low-energy modes. • Ansatz function conditions enforce accurate boundary conditions and convergence. • Consistent linearization enables efficient implicit peridynamic correspondence codes.

Abstract The Peridynamic Petrov–Galerkin (PPG) method is a meshfree approach based on the peridynamic integro-differential form of the momentum equation. The spurious oscillations in the common peridynamic correspondence formulation are investigated. They occur due to an inadmissible linearized mapping of the family deformation field. This leads to a generalized correspondence formulation, which contains the common formulation as a special case. It is based on the weak form of the peridynamic momentum equation. Test and trial function requirements are examined which ensure an exact imposition of Dirichlet and Neumann boundary conditions and Weighted Least Square (WLS) shape functions as well as Local Maximum Entropy (LME) approximants are utilized to examine the PPG Method. A consistent linearization is provided, which can also be used to speed up common implicit peridynamic correspondence codes. It is used in an implicit quasistatic framework to investigate the impact of different shape function combinations. Test cases show that low-energy modes can be prevented by the PPG Method and highlight the fast convergence and stability. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Meshfree method; Petrov–Galerkin method; Peridynamic correspondence formulation; Consistent linearization; Peridynamic reduction; Smoothed particle hydrodynamics

1. Introduction One advantage of meshfree methods is their flexibility in dealing with changing surfaces, large deformations and phase changes. Peridynamics is a non-local integro-differential reformulation of the balance of linear momentum (see Silling [1]). It is widely used in solid mechanics for fracture simulation (see e.g. Silling [2]), but in past years the area of applications rapidly increased (for an overview see Javili et al. [3]). ∗ Corresponding author.

E-mail address: [email protected] (T. Bode). https://doi.org/10.1016/j.cma.2019.112636 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

In the context of Peridynamics two fundamental material modeling approaches are considered as presented in Silling et al. [4]: The first one derives the constitutive response from a non-locally distributed strain energy density function while the second one uses a locally defined strain energy density provided by the classical continuum mechanics. For this purpose, a correspondence formulation and the peridynamic reduction operation provide the link between local and non-local measures. The peridynamic pairwise force state and deformation state are the non-local counterparts to the stress tensor and deformation gradient in the local continuum mechanics theory. With this linking, the resulting non-ordinary state based peridynamic method [5] is a flexible meshfree method that can handle the whole range of local material models. Bessa et al. [6] showed the similarity of the peridynamic correspondence theory to the Reproducing Kernel Particle Method (RKPM) of Liu et al. [7] and other meshfree methods with polynomial basis. In Ganzenm¨uller et al. [8], the peridynamic correspondence formulation was found to be equivalent to the corrected Smoothed Particle Hydrodynamics (SPH) in the total Lagrangian formulation and thus suffers from the same drawbacks. SPH, introduced by Lucy [9] and Gingold and Monaghan [10], is well-established in simulating fluid flows by overcoming the rank deficiency by stabilization techniques, as for example in Monaghan [11]. In terms of peridynamics, the oscillations occur due to the locally averaged deformation gradient, which is a reduction of the peridynamic deformation state (see e.g. Silling [12] and Foster [13]). In recent years, several stabilization techniques were introduced which shall overcome this issue. One approach is, to apply correction forces, see e.g. Littlewood [14],Breitenfeld et al. [15] and Silling [16], which however, leads to an unphysical stiffness of the material. A second technique is based on the stabilization of the displacement field like done in Wu and Ren [17] and Yaghoobi and Chorzepa [18]. Very recently, two concepts were published in which the suppression of low-energy modes is achieved without unphysical stabilization techniques. In Gu et al. [19] the Peridynamic Differential Operator introduced in Madenci et al. [20] is used to calculate the second order derivatives of the displacement field. The imposition of boundary conditions is enforced by means of Lagrangian multipliers. However, this approach is not flexible enough in case of nonlinear material models. The modified correspondence principle of Chowdhury et al. [21] eliminates the lowenergy modes by dividing the family into subdivisions and computes the deformation gradient and resulting force states separately for each subfamily. This method effectively suppresses low-energy modes, whereby an additional treatment in the partitioning into subdivisions is needed. The PPG Method which is developed in the following is aimed to be a completely oscillation-free meshfree method utilizing material models from the local theory without the use of artificial stiffness. The basis for this lies in a comparison between non-local material models based on the deformation state and classical models using the localized deformation gradient. The deformation vector state is a nonlocal measure which contains the deformation of each bond (i.e. the relative position of family particles) and hence is non-linear inside the family H0k . Material models based on this non-linear mapping generally do not suffer from oscillations. As a vector-state generates a more general mapping operation than a linear mapping by a second order tensor, the reduction of the deformation state into a peridynamic deformation gradient yields a loss of information. The assumption that the low-energy modes in the peridynamic correspondence formulation are caused due to this inadmissible linearization leads to the concept of a variable deformation gradient inside a family as demonstrated in Fig. 1 to provide an accurate mapping between reference and current configuration. Thus, for each bond inside the family H0k , including particle k itself, a deformation gradient is computed, respectively. This is analogous to the treatment of underintegration where more integration points include more information. A flexible and fundamental way to achieve this concept consists of a Petrov–Galerkin method that is based on the integro-differential peridynamic linear momentum equation. Different function spaces for the actual and virtual displacement are used to meet their specific requirements. Therefore, a more general correspondence formulation is needed to allow variable deformation gradients inside a family. Furthermore, this formulation paves the way to consistently linearize the residual and to provide an analytical stiffness matrix. This can also be applied within the common peridynamic correspondence models where it is usually computed via Finite Differences which is not efficient. The structure of this paper is as follows: Chapter 2 starts with an introduction to the peridynamic theory and the commonly used correspondence formulation. Subsequent, the derivation of the PPG Method including a generalized correspondence formulation is presented. Shape function requirements and suitable types are presented in chapter 3. A quasistatic framework is derived in chapter 4 which includes the derivation of an analytically computed stiffness

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

3

Fig. 1. The linearized mapping of the deformation of H0k using a constant locally averaged deformation gradient F¯ which is generally non-linear leads to an inadmissible linearized current configuration.

matrix. In chapter 5 an eigenmode analysis proves the absence of low-energy modes for suitable shape functions in contrast to the common peridynamic correspondence formulation. The accuracy and convergence rate of the PPG Method is compared with the Finite Element Method (FEM). Furthermore, the robustness and flexibility is demonstrated in test cases. 2. Peridynamic theory The peridynamic theory is built on the balance of linear momentum. In classical continuum mechanics, it states in a differential form with respect to the reference configuration locally at particle k ρ0k u¨ k = Divk P + ρ0k bk

(1)

with the initial density ρ0k , acceleration u¨ k , first Piola–Kirchhoff stress P and external volumetric forces ρ0k bk . The peridynamic counterpart is written in a non-local integro-differential form in which no spatial derivatives appear as the stress divergence is substituted by an integral over the pairwise interactions within the family H0k : ∫ ( kj ) k k ρ0 u¨ = t − t jk d H0k + ρ0k bk (2) H0k

or in a discretized form using collocation k

ρ0k u¨ k

=

N ∑ (

) tk j − t jk V j cvk j + ρ0k bk ,

(3)

j=1

where the master particle k interacts via the pairwise force densities tk j and t jk with the surrounding neighbor particles j inside the horizon δ k which is set to 2.01 times the particle spacing l¯k in this paper. Generally, a minimal number of neighbor particles is specified by the number of regression coefficients, while large horizons lead to an increasing density of the tangent stiffness matrix and an increased computational time. Here, the sum over all particles inside the family H0k starts with j = 1 which means, that the center particle is excluded. Note that if the center particle is included, the sum starts with j = 0. The neighbor particle volume V j can be multiplied kj by a volume correction factor cv dependent on the particles distance from the family surface of particle k like

4

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

Fig. 2. The geometry of a family is given by the bond vectors ∆Xki and the particle volumes V l inside the horizon δ k . Internal forces acting on particle k due to particle j are considered as the pairwise force densities tk j and −t jk .

e.g. proposed in Madenci and Oterkus [22] to have a more precise family discretization. In this work, the volume correction factor is linearly distributed between zero and one via the following expression: ⎧ ⏐ ⏐ ⏐∆Xk j ⏐ ≤ δ k − l¯k ⎪ 1 ⎪ 2 ⎨ k l¯k ⏐⏐ k j ⏐⏐ ⏐ ⏐ kj δ + 2 −⏐∆X ⏐ k k ¯ ¯ cv = (4) δ k − l2 < ⏐∆Xk j ⏐ ≤ δ k + l2 ⎪ l¯k ⎪ ⏐ ⏐ ⎩ l¯k k j k 0 δ + 2 < ⏐∆X ⏐ . with the position difference vectors (called bonds) in the reference configuration ∆Xk j . These are usually stored in a container, a so-called vector state, as follows: ⎛ ⎞ ⎛ ⎞ ∆Xk1 ∆xk1 ⎜ ∆Xk2 ⎟ ⎜ ∆xk2 ⎟ ⎜ ⎟ ⎜ ⎟ ∆Xk j = X j − Xk → ∆Xk = ⎜ or ∆xk j = x j − xk → ∆xk = ⎜ (5) .. .. ⎟ ⎟. ⎝ ⎠ ⎝ ⎠ . . ∆Xk N

k

∆xk N

k

The bond vector state of particle k in the current configuration is described by xk . An illustration of the discretized neighborhood of particle k with the kinematic and dynamic connections to its family H0k is provided in Fig. 2. The principal of pairwise force densities enables the conservation of linear momentum in the discretized form between two particles and independently on the pairwise force expression: ( ) ( ) ∂Ik j ∂I jk = tk j − t jk V k V j cvk j = − t jk − tk j V j V k cvjk = − . (6) ∂t ∂t where it is provided that δ k = δ j and l¯k = l¯ j . For a detailed introduction in the theory of states and its operations the reader is referred to Silling et al. [4]. The linking to the first Piola–Kirchhoff stress tensor from local continuum mechanics theory is now performed via a correspondence formulation, which can be derived via a comparison of local and non-local strain energy densities like in Silling et al. [4] or Madenci and Oterkus [22]. 2.1. Correspondence formulation The commonly used correspondence formulation, first introduced in Silling et al. [4], states tk j = ωk j Pk · Kk

−1

· ∆Xk j .

The locality of the correspondence formulation can be adjusted with a radial weight function ω on the length of the reference bonds. In the current work, a radial weight function (kernel) of 1 ωk j = ⏐ ⏐ ⏐∆Xk j ⏐2

(7) kj

which depends (8)

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

5

seemed to be an appropriate choice for most stable results. The symmetric shape tensor Kk (also called mass matrix or correction matrix) occurring in the correspondence formulation of Eq. (7) is a measure for the shape of a family in the reference configuration and is defined as the tensor product of the reference bond vector states: k

k

k

k

∫

ω ∆X ⊗ ∆X kj

K = ∆X ∗ ∆X = H0k

kj

kj

d H0k

=

N ∑

ωk j ∆Xk j ⊗ ∆Xk j V j cvk j .

(9)

j=1 k

The first Piola–Kirchhoff stress Pk at particle k depends on the locally averaged deformation gradient F , which can be obtained by a reduction of the current bond vector state: ⎛ k ⎞ ⎛ k ⎞−1 N N ∑ ∑ ( ) k −1 F = R{∆xk } = ∆xk ∗ ∆Xk · Kk = ⎝ ωk j ∆xk j ⊗ ∆Xk j V j cvk j ⎠ · ⎝ ωk j ∆Xk j ⊗ ∆Xk j V j cvk j ⎠ j=1

j=1

(10) This can now be used for a constitutive equation like the Neo-Hookean material (see e.g. Wriggers [23]) to compute a family-averaged stress: ( ) ( ) −1 −1 Sk = Λk J k J k − 1 Ck + µk 1 − Ck , (11) where Sk is the second Piola–Kirchhoff stress, Λk and µk are material parameters, J k is the determinant of the locally averaged deformation gradient and kT

Ck = F

k

·F

(12)

is the right Cauchy–Green tensor. A conversion to the first Piola–Kirchhoff stress closes the cycle to the correspondence formulation equation (7): k

Pk = F · Sk .

(13)

2.2. Peridynamic Petrov–Galerkin method: A generalized correspondence formulation In the correspondence formulation of Eq. (7), the pairwise forces depend only on the stress at particle k and thus only on the locally averaged deformation gradient at particle k. As already mentioned, this correlates with a loss of information towards the peridynamic deformation vector state. Hence, to enable the benefit of a variable deformation gradient field inside a family, a more general correspondence formulation is needed which includes changing stresses. In the following, the index notation used in the PPG Method is presented. The expression Fk j will stand for the deformation gradient at the position of particle j referring to the family of particle k as the deformation gradient within a family is non-constant (see Fig. 1). Note, that Fk j is generally not the same as F j j , even if it denotes the deformation gradient at the position of particle j the reference differs. In contrast, in the case of classical Peridynamics Fk j is equal to Fkk as there is a constant deformation gradient per family. Furthermore, family independent quantities like the volume V k are written with a single superscript which indicates the concerning particle. Differences like the reference position difference between particle k and j are denoted like before as ∆Xk j = X j − Xk . In case of shape functions, a triple index is needed. Therefore N k ji stands for the ith shape function of the family of particle k evaluated at the position of j. In the following, the generalized correspondence formulation is derived by comparing the weak forms of the local and non-local conservation of linear momentum equations (1) and (3). This is similar to the derivation of the common peridynamic correspondence formulation of Madenci and Oterkus [22]. First, the peridynamic form has to be multiplied by a virtual displacement δuk and integrated over the entire domain Ω0 to obtain the weak form: ∫ ∫ ∫ ∫ ( kj ) k k k jk k k ρ0 u¨ · δu dΩ0 = t − t · δu d H0 dΩ0 + ρ0k bk · δuk dΩ0 . (14) Ω0

Ω0

H0k

Ω0

6

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

The weak form of the local momentum equation (1) (see for instance Wriggers [23]) yields ∫ ∫ ∫ ∫ ρ0k u¨ k · δuk dΩ0 = − Pk : δFk dΩ0 + Tk · δuk d∂Ω0 + ρ0k bk · δuk dΩ0 . Ω0

∂Ω0

Ω0

(15)

Ω0

Considering the fact, that the pairwise force densities tk j act only within the family H0k of particle k, the integral over d H0k can be extended over the whole body dΩ0′ . Hence the strain energy term of the peridynamic weak form (14) can be rewritten as ∫ ∫ ∫ ∫ ∫ ∫ ( kj ) t − t jk · δuk d H0k dΩ0 = tk j · δuk dΩ0′ dΩ0 − t jk · δuk dΩ0′ dΩ0 k ′ ′ Ω0 H0 Ω Ω Ω Ω ∫ 0∫ 0 ∫ 0∫ 0 jk j ′ = t · δu dΩ0 dΩ0 − t jk · δuk dΩ0′ dΩ0 ′ ′ Ω Ω Ω0 Ω0 ∫ 0∫ 0 (16) jk kj k = t · ∆δu d H0 dΩ0 Ω0 H0k ∫ =− Pk : δFk dΩ0 , Ω0

where the indices k and j are changed in the second line and the integration domain is reset to d H0k in the third. By using the following ansatz for the virtual displacements k

N ( ) ∑ δu ∆Xkk = Nηkk j ∆δuk j ,

(17)

j=0 k ji

where Nη

δFkk

is the ith test function of H0k at the position of j, the virtual deformation gradient follows ( ( )) Nk kk j ∑ ∂u ∆Xkk ∂ Nη = ∆δuk j ⊗ . =δ 1+ ∂X ∂X j=0

(18)

Equating the strain energy in the peridynamic and local form (16) and applying the virtual deformation gradient (18) yields ∫

k

N ∑

t

jk

kj

· ∆δu V

Ω0 j=1

j kj cv

∫ dΩ0 = −

k

N ∑

Ω0 j=1

1

kk j

Pk · kj

V j cv

∂ Nη · ∆δuk j V j cvk j dΩ0 . ∂X

(19)

Thereby, the generalized correspondence formulation follows kk j

1

jk

∂ Nη · ∂X

j jk

∂ Nη t =− P or rather t =− P · . (20) kj jk j k ∂X V cv V cv Since the deformation gradient is now non-constant within a family, the approximation at an absolute position depends on the origin of the family. Thus, also the first Piola–Kirchhoff stress depends on this reference and this is why the pairwise force density and its stress counterpart have to have the same origin. Inserting the correspondence formulation (20) back into the strong form of the peridynamic momentum equation (3) yields the momentum equation of the PPG Method: ( ) Nk kk j kj j jk ∑ V j cv k j ∂ N η jk ∂ Nη k k P · (21) ρ0 u¨ = − P · + ρ0k bk . k c jk ∂X ∂X V v j=1 jk

kj

1

kj

Furthermore, an ansatz for the primary variables is needed to compute the deformation gradient which is needed for the local constitutive equation like for example the Neo-Hooke material (see Eq. (11)): k

N ( ) ∑ u ∆Xk j = Nuk ji ∆uki , i=0

(22)

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

7

k ji

with the trial functions Nu of H0k . By computing the gradient and adding the second order unit tensor 1, the deformation gradient at particle j with respect to the family H0k follows k

N ∑

Nk

k ji

k ji

∑ ∂ Nu ∂ Nu F =1+ ∆u ⊗ = . ∆xki ⊗ ∂X ∂X i=0 i=0 kj

ki

(23)

The above Eqs. (21) and (23) are quite general. The peridynamic correspondence formulation (7) and locally averaged deformation gradient (10) can be restored by using specific test and trial functions as demonstrated in Appendix A. 2.3. Neumann boundary conditions By investigating the classical and peridynamic weak forms (15) and (14) it stands out that the stress boundaries on the surface do not have a counterpart in the peridynamic form. This is why, external forces are converted into volumetric forces in the first particle layer. In this way, the external volumetric forces consist of two parts, truly volumetric ones like the gravity and converted ones: ∫ ∫ ∫ qk ρ0k gk · δuk dΩ0 + ρ0k bk · δuk dΩ0 = · δuk dΩ0 , (24) k A Ω0 Ω0 Ω0 where gk is the gravity, qk is a line load acting at particle k and Ak is the surface of particle k. 3. Shape functions The following chapter deals with the choice of suitable test and trial functions. First, requirements that have to be imposed by the different shape functions are presented. Thereupon, two classes of shape functions are introduced. Referring to Krongauz and Belytschko [24], a key ingredient for convergence is consistency of the trial functions. Linear consistency for instance, guarantees the exact approximation of a linear function. Mathematically, this consists of the following conditions on the trial function: The zeroth order consistency condition with its spatial derivative states k

N ∑

k

Nuk ji

= 1,

i=0

N k ji ∑ ∂ Nu i=0

∂X

= 0.

(25)

First order consistency holds if k

N ∑

k

Nuk ji ∆Xki = ∆Xk j ,

i=0

N ∑

k ji

∆Xki ⊗

i=0

∂ Nu =1 ∂X

(26)

and second order consistency is achieved analogously: k

N ∑

Nuk ji ∆Xki ⊗ ∆Xki = ∆Xk j ⊗ ∆Xk j .

(27)

i=0

Higher consistency is supposed to lead to better convergence rates. In a non local method like the PPG Method, a linear regression polynomial is not sufficient to prevent low energy oscillations. The subsequently listed weighted least square shape functions are consistent up to the order of the regressed polynomial, and thus a good choice for meshless trial functions. However, the requirements for the test functions are different. According to the principle of virtual displacements, only geometrically possible virtual displacements are allowed. Therefore, the virtual displacement has to be zero at fixed displacement boundaries. One possibility to impose these Dirichlet boundaries is the fulfillment of Kronecker-δ property at the particles in question: N k ji = δ ji .

(28)

8

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

3.1. Peridynamic integration constraint Gauss’s theorem leads to a further condition for accurate modeling of free and bounded surfaces. Applying it on the inner forces, the stress divergence from Eq. (1) integrated over particle k yields ∫ ∫ k k Div P d V = P · Nk d A k . (29) Vk

Ak

Exchanging the left side by its peridynamic counterpart (21) and setting the stress to be constant leads to the peridynamic integration constraint { Nk j jk kj kk j ∑ 0 on Ω0 j cv ∂ N η k k k ∂ Nη − V jk =N A = , (30) V k k A N on ∂Ω0 ∂X cv ∂X j=0 where Ak is the surface of particle k which only exists, if the particle is at the surface of the entire body Ω0 . The peridynamic integration constraint ensures the accurate imposition of Neumann boundaries as described in subchapter 2.3. Note, that the left term of the sum over j is equal to zero if the test functions are zeroth order consistent (see Eq. (25)). Further requirements for the test function consist in existent derivatives at the given data points (i.e. family particles) and a good numerical integrability. 3.2. Weighted least square shape functions Motivated by the method of least squares (see for instance Lancaster and Salkauskas [25]), this subsection deals with the preparation of polynomial regression shape functions that fit into the shape function approach of Eq. (22) which is based on differences. Consider given scalar values f j at the neighbor particles and a value f k at the master particle. The difference ∆ f k j is defined as ∆ f kj = f j − f k .

(31)

The regression ∆ fˆk of the scalar state ∆ f k inside H0k is gained from Taylor’s expansion neglecting the constant term as we only deal with differences: ∆ fˆk (∆X) = f k j = a1k ∆X + a2k ∆Y + a3k ∆X 2 + a4k ∆X ∆Y + a5k ∆Y 2 + · · · = ak · p ,

(32) k

where the monomials are stored in a vector p and the currently unknown coefficients in a vector a . Note, that the monomials do not generally satisfy point symmetry as in the linear case (i.e. pk j ̸= −p jk ). Since ∆ fˆk is a regression function, an error exists at the data points: ( ) ϵ k j := ∆ fˆk ∆Xk j − ∆ f k j . (33) Using least square fitting, the best fitted function can be gained by minimizing the sum of weighted error squares. Here, a radial weight ωki (also called kernel or window function) and a volume weight V i is applied: k

E k :=

N ∑

( )2 ωki ϵ ki V i cvki .

(34)

i=0

The minimization problem can also be stated as an overdetermined linear system of equations in state notation ∆ f k = pk · ak .

(35)

By tensor multiplying the shape vector state pk and multiplying the inverse of the resulting shape tensor Mk from the left one obtains the coefficients ak as Nk ( ) ∑ k k k −1 k k −1 a =M · p ∗∆f =M · ωki pki ∆ f ki V i cvki , (36) i=0

with the symmetric higher order shape tensor (also called correction or mass matrix) k

Mk = pk ∗ pk =

N ∑ i=0

ωki pki ⊗ pki V i cvki .

(37)

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

9

Applying the now known coefficients ak from Eq. (36) in the regression function (32) yields k

k

N N ∑ ( ) ∑ ( ) −1 ∆ fˆk ∆Xk j = ωki V i cvki pki · Mk · p ∆Xk j ∆ f ki = N k ji ∆ f ki , i=0

(38)

i=0

where the regression function ∆ fˆk is rewritten in terms of ansatz functions. N k ji stands for the ith shape function of the family H0k , evaluated at the reference position ∆Xk j of particle j. Here, the peridynamic Weighted Least Square (WLS) shape functions for 1 ≤ i ≤ N k follow N k ji = ωki V i cvki pki · Mk

−1

· pk j .

(39)

By evaluating the derivatives of the shape functions at the data points, the computation simplifies to: ( ) ∂pk j ∂ N k ji −1 = ωki V i cvki pki · Mk · . (40) ∂X ∂X As the left and right terms of Eq. (38) are zero for i = 0, the shape function follows via the partition of unity which is included by the linear consistency conditions. Hence for i = 0, the peridynamic regression shape function and its derivative are k

k

N

k jk

=1−

N ∑

N

k ji

N ∑ ∂ N k ji ∂ N k jk =− . ∂X ∂X i=1

and

i=1

(41)

Considering a linear displacement field, each regression polynomial which has a linear term is able to approximate the field exactly, i.e. the sum of the weighted error squares E k = 0 and all coefficients of higher order terms go to zero. Therefore, linear consistency is enforced by polynomial regression shape functions of first order. Analogously, a shape function of nth order can exactly approximate an nth order displacement field. Remark. Analogously to the above derivation of regression and its derivative of a scalar field, a vector and even higher valued regression can be defined, whereby the order of the coefficient matrix increases. When applied to the peridynamic theory, the derivative of this regression function can be called a peridynamic reduction. Hence, the nth order reduction of the vector state ∆fk yields ( ) −1 ∂p (∆X) n R{∆fk } (∆X) = ∆fk ∗ pk · Mk · . (42) ∂X Note, that the locally averaged deformation gradient from Eq. (10) which is used in the peridynamic correspondence formulation can be gained by the use of linear regression where p = ∆X: ( ) ( ( ) ) k −1 ∂p (0) −1 ∂∆X −1 1 = ∆xk ∗ ∆Xk · Kk · = ∆xk ∗ ∆Xk · Kk = F . (43) R{∆xk } (0) = ∆xk ∗ pk · Mk · ∂X ∂X 3.3. Local maximum entropy shape functions Besides the polynomial regression shape functions, there exist several different types of shape functions for an arbitrary number of data points. One of these is the Local Maximum Entropy (LME) approximation introduced by Arroyo and Ortiz [26] which possesses a weak Kronecker-δ property and provide linear consistency. A detailed investigation of LME shape functions including a two-dimensional visualization can be found in Weißenfels and Wriggers [27]. LME approximants are constructed from normalizing radially exponential functions fulfilling the first order consistency condition by means of Lagrangian multipliers: ⏐ ⏐2 ( ) −β k ⏐∆Xk j −∆Xki ⏐ +λk Z ki ∆Xk j = e ⏐

with the normalization ( ) Z ki ∆Xk j k ji N =∑ k ( ). N ki ∆Xk j Z i=0

⏐

(

)( ) ∆Xk j · ∆Xk j −∆Xki

(44)

(45)

10

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

The parameter β k controls the locality of the shape functions. In order to gain independence from the dimensions it is stated as βk =

γk l

(46)

k2 k

where the parameter γ k can be chosen between 50 and 100 for horizon sizes of δ k ≥ 1.5l and the characteristic k length l is computed as twice the horizon: k

l = 2δ k .

(47)

The first order consistency conditions is ensured by enforcing the Lagrangian multipliers λk to hold k

r ∆X , λ k

(

kj

k

)

=

N ∑

N k ji ∆X ji = 0 .

(48)

i=0

The Newton–Raphson algorithm (see e.g. Wriggers [23]) can be used to solve the non-linear equation (48) and the shape function derivatives can be obtained via an automatic differentiation tool as for instance AceGen (see also Korelc and Wriggers [28]). Even if the LME shape functions do not fulfill the test function criterion (30) exactly and the required Kronecker-δ property is only met at the convex boundary, these are closely fulfilled for strongly localized LME shape functions and therefore lead to good results. 4. Implicit quasistatic framework In this chapter, the framework in which the PPG Method is validated and tested is presented. In a quasistatic simulation the external forces are applied in load steps. Dirichlet boundary conditions are treated in terms of an additional layer of displacement constraint particles that contribute to the pairwise force densities and therefore implement the reaction forces. In each load step the Newton–Raphson algorithm is used to find the equilibrium position in which internal and external forces cancel each other. Considering constant quantities inside each particle, the peridynamic weak form (14) can be written as { ( )} ∫ ( kj ) ne k k jk k k k k t −t d H0 − ρ0 V b = 0. (49) G (u, δu) = ∪k=1 δu · −V H0k

In each Newton–Raphson iteration, the residual is linearized with respect to displacements which are then corrected by solving the formed linear system of equations. After each iteration, the residual is recomputed and checked for convergence. In this way, one iteration step consists of two main parts. On one hand, there is the computation of a tangent matrix KΩ0 where its accuracy is crucial for achieving a good convergence rate in the Newton algorithm. In the field of Correspondence Peridynamics, this tangent stiffness matrix is usually approximated as a secant via Finite Differences, which is costly and results, according to Brothers et al. [29], in the majority of the computation time in implicit peridynamic analyses. The PPG Method provides a much faster way to assemble the tangent by consistently linearizing the residual even for higher order approaches where the constitutive law has to be evaluated at every neighbor particle. On the other hand, the linear system of equations has to be solved. Note, that the tangent matrix is generally not symmetric. In this work, the biconjugate gradient stabilized method (BiCGSTAB) is used with an Incomplete LU preconditioner. 4.1. Secant matrix via finite differences Utilizing the Central Finite Differences to obtain a stiffness matrix is a relatively simple way to use the broad spectrum of materials from the local theory in implicit peridynamic frameworks (see e.g. Littlewood [30]). In doing so, each particle has to be perturbed in every dimension in positive and negative direction. The change in the residual of influenced particles has to be computed and in the global stiffness matrix. Note, that the ( assembled ) k kj influenced particles are located in a sphere of radius 2 δ k + l2 or rather 2δ k if no volume correction cv is used.

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

11

Table 1 Petrov–Galerkin shape function combinations. Combination

Trial function

Test function

L-L L-LME pL-pL Q-L Q-LME

First order WLS First order WLS Piecewise first order WLS Second order WLS Second order WLS

First order WLS LME approximants Piecewise first order WLS First order WLS LME approximants

The concerning entry of the residual change of particle k in x-direction while perturbing particle i in y-direction states ( ) ( ) r xk u + ϵ yi − r xk u − ϵ yi Ω0 K (k + x, i + y) ≈ . (50) 2ϵ The perturbation size ϵ should be chosen as small as possible to get the best approximation, but not too small such that numerical errors dominate. Littlewood [30] suggests a value of about 10−6 times the particle spacing. 4.2. Tangent matrix via linearization In this paragraph, the tangent stiffness matrix is derived via consistent linearization for materials from the local theory and is used for the test cases, see chapter 5. For complex materials automatic differentiation tools as for instance AceGen (see also Korelc and Wriggers [28]) can be employed. Considering only displacement-independent loads and a non changing discretization, the residual of Eq. (49) is linearized as follows: } { ∫ ∫ ( kj ) ne jk k k k . (51) t − t · δu d H0 d V ∆G = ∪k=1 ∆ − Vk

H0k

In the further treatment of ∆G, it can be split into pairwise geometric and substantial parts. The detailed linearization can be found in Appendix B. The global tangent stiffness matrix can be assembled as follows: ⎫ ⎧ Nk ∑ Nk ( Nj ( Nk ∑ ⎨∑ )⎬ ) ∑ ij ij ij ij e e G 2 1 + M2 . (52) Kk = ∪nk=1 G 1 1 + M1 + KΩ0 = ∪nk=1 ⎭ ⎩ j=1 i=0

j=1 i=0

Note that the inner sum of the right pairwise part operates on the family of particle j from the outer sum. This leads to a rectangular particle-wise stiffness matrix Kk . 5. Numerical examples The methodology developed in the previous chapters is now tested and investigated with the implicit quasistatic PPG framework. First, an eigen mode analysis is performed to prove the absence of low-energy modes for suitable shape functions. Second, Cook’s membrane is used to compare different shape function types (see Table 1) regarding the existence of spurious oscillations and the accurate imposition of displacement boundaries. Further, the convergence rates are compared to FEM by means of a punch test. The robustness even for large deformations is exemplarily demonstrated. Finally, the flexibility of the PPG Method is shown in a contact simulation. 5.1. Eigen modes analysis As the particle-wise stiffness matrices are rectangular, it is not possible to compute eigenvalues of a family. Therefore, the global tangent stiffness matrix of chapter 4.2 is used to get the eigen values and modes of a bloc of rubber. Independent on the choice of shape functions, there exist three zero eigen values in 2-D — two translational and one rotatory. In case of linear weighted least square shape functions, as in the common peridynamic correspondence theory, there exist a bunch of spurious low-energy modes before the first normal eigen mode as demonstrated in Fig. 3(a). By simply making a quadratic regression of the deformation field inside a family, i.e. using second order weighted least square shape functions as trial functions, the low-energy modes disappear (see Fig. 3(b)) and the physical eigen modes are in good agreement to those obtained with FEM.

12

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

Fig. 3. Non-zero eigen modes of a rubber bloc using linear or higher order weighted least square trial functions or Q2 finite Elements.

Fig. 4. Geometry, material data and boundary conditions of Cook’s membrane benchmark.

5.2. Cook’s membrane: Investigation of spurious oscillations and Dirichlet boundaries The behavior of different shape function combinations of Table 1 for a nearly incompressible Neo-Hookean material is investigated with Cook’s Membrane which is a standard benchmark for Finite Elements (see for example Kadapa et al. [31]). The geometry, boundary conditions and material data are visualized in Fig. 4. In a quasistatic simulation, the load q0 is applied linearly in 2 load steps for each combination, respectively.

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

13

N Fig. 5. Hydrostatic pressure fields after applying the load of q0 = 6.25 m in 2 load steps in a quasistatic simulation using the different shape function types from Table 1.

Fig. 6. On the left: Spurious large deformations near the boundary due to fitting test functions. In the middle: Smooth transition from wall particles (red) to surrounding free particles (blue) by account of weakly Kronecker-δ preserving test functions. On the right: Finite element solution with Q2 elements. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The hydrostatic pressure is plotted in Fig. 5. It can be seen, that in the case of linear regression test and trial functions (i.e. common Peridynamics), low-energy modes dominate the displacement and pressure field. A quadratic weighted least square approach for the trial function (Q-L and Q-LME) leads to very smooth displacement and pressure fields without any oscillations. Subsequently, the load is increased to q0 = 50 mN to reveal the inaccuracy of Kronecker-δ violating test functions regarding Dirichlet boundaries. In Fig. 6, there is the shear stress colored and a zoom of the upper left region of Cook’s membrane. While the linear regression test functions cause non-physical large displacements near the wall particles, the Q-LME combination satisfies a smooth transition from fixed to free particles due to the more precise enforcement of Dirichlet boundaries by means of weakly Kronecker-δ preserving shape functions. The stress distribution shows good agreement to the FEM Q2 reference solution.

14

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

Fig. 7. On the left: Geometry, material data and boundary conditions of the punch test. The wall particles on the right and bottom are fixed in both coordinate directions whereas on the left they are just horizontally fixed. On the right: Exemplary end configuration for the Q-LME combination with a particle spacing of 2.5 mm.

Fig. 8. Convergence analysis: The minimal vertical displacement is plotted over the number of particles. Different PPG combinations from Table 1 are compared with FEM solutions of bi-linear (Q1) and bi-quadratic quadrilaterals (Q2).

5.3. Punch test: A convergence analysis After finding that second order consistent shape functions are suitable to suppress low-energy modes, the convergence behavior is now compared to FEM using a punch test. Here, a block of rubber is investigated which is compressed as depicted in Fig. 7. The load is applied in 2 load steps in a quasistatic framework and the shape function combinations are compared with varying spatial discretizations. The final state of the rubber block is exemplarily shown in Fig. 7 for the Q-LME combination. The convergence analysis is performed by plotting the minimal occurring vertical displacement over the number of nodes (see Fig. 8). As one can see, the PPG combinations using polynomial regression test functions converge

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

15

Fig. 9. Stable high deformation punch test using a Q-LME with a nodal spacing of ∆l = 0.5 mm and a horizon of δ = 2.01∆l.

relatively bad due to the non-accurate enforcement of Dirichlet boundary conditions. The LME approximants fulfill the test function conditions of Section 3 much better. Hence, it converges faster to the FEM solution where the convergence rate is between those of FEM with Q1 and Q2. 5.4. Punch test: Large deformation stability Subsequently, the most exact and fast converging Q-LME combination is used to demonstrate the stability of the PPG Method even for very large deformations. Here, the punch test from the previous subsection is amended by increasing the load by a factor of 30. To prevent problems due to non-linear force boundaries, the horizontal positions of the concerning particles are fixed. The end configuration of the quasistatic loading is visualized in Fig. 9, where the vertical displacement is colored. The deformation is smooth and completely free of spurious oscillations. 5.5. Punch test with penalty contact The last numerical example includes a punch test where a block of rubber is compressed by a barrier as illustrated in Fig. 10. Referring to Ganzenm¨uller [32], these problems are hard to simulate using FEM. In this paper, the constraint of the barrier is enforced using the penalty method. The barrier is incrementally shifted from u = 0 mm up to u = 0.4 mm in a quasistatic simulation using the Q-LME combination of the PPG Method. The resulting configurations in case of four different Poisson’s ratios are shown in Fig. 11, where the shear stress is colored. At the right end of the barrier, the stress is heavily localized but has no unphysical patterns. 6. Conclusion In this work the spurious low-energy modes occurring in the peridynamic theory of correspondence materials are found out to stem from an inadmissible linearized family-mapping of the reference configuration to the current configuration by a locally averaged deformation gradient. The PPG Method which is presented paves the way for varying deformation gradients inside a particles family and prevents the presence of low-energy modes. For future research, the concept of non-linear mapping inside a family is extendible to a diverse field of applications.

16

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

Fig. 10. Geometry, material data and boundary conditions of the punch test. The wall particles on the left are horizontally and at the bottom vertically fixed. The upper particles interact with a penalty barrier if their x-position is smaller than 3 mm.

Fig. 11. End position of barrier contact punch test with varying Poison’s ratios. The shear stress is strongly localized at the end of the penalty barrier.

The investigation of test and trial function requirements provides an explanation for the accuracy of the shape functions used in the numerical test cases. The quadratic WLS shape function is suitable as a trial function due to its second order consistency. On the other hand, the LME approximants match the test function conditions in an optimal way which is reflected in the numerical test. In chapter 4 the quasistatic PPG framework is presented which is used in the numerical studies. Besides a Finite Difference approximation of the tangent stiffness matrix, a much faster way of assembling the stiffness matrix via a consistent linearization of the residual is presented. This procedure can be used to improve common peridynamic correspondence codes. Subsequently, from several test cases the PPG Method was found to be a stable and fast converging meshfree method which does not need stabilization techniques, provided that suitable shape function combinations are used. Note, that the common Peridynamics with correspondence materials and thus also the total Lagrangian Smoothed Particle Hydrodynamics is actually a special case of linear regression shape functions (see Appendix A). Appendix A. The common correspondence formulation as a special Case of the PPG method The equivalence of the common peridynamic correspondence formulation with the PPG Method can be obtained by the following test and trial function derivatives: j jk

∂ Nη ∂X

= ωk j V k cvjk Kk

−1

k ji

· ∆X jk

and

∂ Nu −1 = ωki V i cvki Kk · ∆Xki . ∂X

(53)

17

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

Inserting these in the generalized correspondence formulation of Eq. (20) yields the classical one of Eq. (7): tk j = − =−

1

j jk

Pk j · jk

V k cv 1

jk V k cv kj kj

∂ Nη ∂X

Pk j · ωk j V k cvjk Kk

= ω P · Kk

−1

−1

(54)

· ∆X jk

· ∆Xk j ,

where the first Piola–Kirchhoff stress is as well as the deformation gradient constant in a family. The deformation gradient states by applying the trial functions of (53) on Eq. (23): k

kj

F =

N ∑

k ji

∆xki ⊗

i=0

∂ Nu ∂X

k

=

N ∑

∆xki ⊗ ωki V i cvki Kk

−1

· ∆Xki

(55)

i=0 k

=

N ∑

ωki ∆xki ⊗ ∆Xki V i cvki · Kk

−1

,

i=0

which is equivalent to the locally averaged deformation gradient of Eq. (10). Appendix B. Derivation of the linearized stiffness matrix Taking into account that the pairwise force densities are defined to be zero outside of H0k , Eq. (51) can be written as ∫

∫

(

∆G = −∆ Ω0

H0k

) tk j − t jk · δuk d H0k dΩ0 = −∆

Changing the indices k and j yields ∫ ∫ ∫ ( jk ) ∆G = −∆ t − tk j · δu j dΩ0′ dΩ0 = ∆ Ω0

Ω0′

∫

∫

Ω0

∫ Ω0

( Ω0′

( H0k

) tk j − t jk · δuk dΩ0′ dΩ0 .

) tk j − t jk · δu j d H0k dΩ0 .

(56)

(57)

It follows the family-wise view under assumption of constant initial particle volumes and actual displacement independent virtual displacements: ∫ ∫ ∫ ∫ ( ) ( kj ) jk j k k k (58) t − t · δu d H0 d V = ∆ tk j − t jk · δu j d H0k d V k . ∆G = ∆ Vk

H0k

Vk

H0k

By inserting the correspondence formulation of Eq. (20) it changes to ( ) ∫ ∫ j jk kk j 1 1 k k j ∂ Nη jk ∂ Nη ∆ − ∆G = P · + P · · δu j d H0k d V k . jk kj k j k k ∂X ∂X V cv V cv V H0

(59)

Applying center particle integration on the outer integral and neighbor particle integration on the inner integral kj jk considering cv = cv it leads to ) ( Nk j jk kk j ∑ k j k j ∂ Nη k jk ∂ Nη ∆G = ∆ −V P · +V P · · δu j . (60) ∂X ∂X j=1 By using the product rule and taking into account that the test functions are defined in the reference configuration it simplifies to ( ) Nk j jk kk j ∑ k j k j ∂ Nη k jk ∂ Nη ∆G = −V ∆P · + V ∆P · · δu j . (61) ∂X ∂X j=1

18

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

Substituting the first Piola–Kirchhoff stress for P = F · S and the linearized second Piola–Kirchhoff stress for ∆S = C : ∆E yields via the product rule ( Nk ∑ ( ) ∂ Nηj jk k ∆G = +V j ∆Fk j · Sk j + Fk j · Ck j : ∆Ek j · ∂X j=1 ) (62) ( ) ∂ Nηkk j k jk jk jk jk jk −V ∆F · S + F · C : ∆E · · δu j ∂X = T1 + T2 + T3 + T4 ( ) with the fourth order material tangent C and the Green–Lagrange strain E = 21 FT · F − 1 . The further breakdown of the above expression is separated into a geometric (T1 and T3 ) and a substantial part (T2 and T4 ). In case of the Neo-Hookean material from Eq. (11), the material tensor states ( 2 )] ( −1 ) ( ) 1[ j −1 −1 −1 −1 −1 kj Cabcd = Λ j J k j 2J k j − 1 C k j ab C k j cd + 2µ − Λ j J k j − 1 C k j ac C k j bd + C k j ad C k j bc . (63) 2 B.1. Geometric part Using the test function approach of Eq. (22) and the derivative of the partition of unity condition, the linearized deformation gradient ∆Fk j gets k

N ∑

k ji

Nk

k ji

∑ ∂ Nu ∂ Nu ui ⊗ = . ∆F = ∆H = ∆u ⊗ ∂X ∂X i=0 i=0 kj

kj

ki

(64)

Thereby, the first geometric part writes k

k

T1 =

N N ∑ ∑

k ji

−V j ui ⊗

j=1 i=0

j jk

∂ Nu ∂ Nη · Sk j · · δu j . ∂X ∂X

Rearranging this leads to ( ) Nk Nk ∑ Nk Nk ∑ k ji j jk ∑ ∑ ij j j ∂ Nu k j ∂ Nη δu · −V δu j · G 1 1 · ui T1 = ·S · 1 · ui = ∂X ∂X j=1 i=0 j=1 i=0

(65)

(66)

with the geometrical scalar k ji

j jk

∂ Nu ∂ Nη · Sk j · . (67) ∂X ∂X Analogously, the second geometric part T3 is treated by changing the indices of the linearized deformation gradient (Eq. (64)): ij

G 1 = −V j

k

T3 =

j

N ∑ N ∑

jki

V k ui ⊗

j=1 i=0

kk j

∂ Nη ∂ Nu · S jk · · δu j . ∂X ∂X

Again, a rearrangement leads to ( ) Nj Nk ∑ Nk ∑ Nj jki kk j ∑ ∑ ∂ N ∂ N η u ij j k jk i δu · V T3 = ·S · 1·u = δu j · G 2 1 · ui ∂X ∂X j=1 i=0 j=1 i=0

(68)

(69)

with jki

ij

G2 = V k

kk j

∂ Nu ∂ Nη · S jk · . ∂X ∂X

(70)

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

19

B.2. Substantial part By using the linearized deformation gradient of Eq. (64) and considering the product rule, the linearized Green–Lagrange strain can be written as ( ) Nk k ji k ji ) ∑ 1 ( kj T kj ∂ Nu 1 ∂ Nu kj kj T kj i kj kj T i ∆E = ∆F · F + F · ∆F = ⊗u ·F +F ·u ⊗ . (71) 2 2 ∂X ∂X i=0 Inserting this into the first substantial part T2 leads to ( ) Nk ∑ Nk k ji j jk ∑ ∂ Nu ∂ Nη 1 j kj kj i kj T2 = ⊗u ·F · · δu j − V F ·C : 2 ∂X ∂X j=1 i=0 ( ) k Nk N k ji j jk ∑∑ ∂ Nu ∂ Nη 1 j kj kj T i kj · · δu j . + − V F ·C : F ·u ⊗ 2 ∂X ∂X j=1 i=0

(72)

Written in index notation and considering the symmetries of C, the two parts can be merged as follows: k

k

N ∑ N ∑ j=1

k ji

Nk

j jk

Nk

k ji

j jk

∑∑ 1 ∂ Nu i k j ∂ Nη ∂ Nu ∂ Nη 1 kj kj j kj kj kj T j u e Fed δu f + δu f − V j F f a Cabcd − V j F f a Cabdc Fde u ie 2 ∂ X ∂ X 2 ∂ X ∂ X c b c b i=0 j=1 i=0 k

=

k

N N ∑ ∑ j=1

j jk

Rewriting this in matrix notation yields [ ( ) ] Nk Nk Nk ∑ Nk ∑ k ji j jk ∑ ∑ ∂ Nη ij j j kj k j ∂ Nu kj T i δu · −V F · ·C · ·F δu j · M1 · ui T2 = ·u = ∂X ∂X j=1 i=0 j=1 i=0 with the substantial matrix ( ) j jk k ji ∂ N ∂ N T η u ij M1 = −V j Fk j · · Ck j · · Fk j . ∂X ∂X Utilizing the Voigt notation for Eq. (75), it can be written as ( ) T T ij M1 = −V j Fk j · Bηj jk · Dk j · Bku ji · Fk j with the B matrix in 2-D ⎛ k ji ⎞ ∂N 0 ∂X ⎜ ∂Nk ji ⎟ . Bk ji = ⎝ 0 ⎠ ∂Y ∂Nk ji ∂Y

(73)

k ji

∂ Nu j k j ∂ Nη kj kj T −δu f V j F f a Cbadc Fde u ie . ∂ Xb ∂ Xc i=0

(74)

(75)

(76)

(77)

∂Nk ji ∂X

Analogously to the first substantial part, the second one can be handled by changing the indices k and j of the linearized Green–Lagrange strain tensor (Eq. (71)) and inserting it into T4 : ( ) Nk ∑ Nj jki kk j ∑ 1 k jk ∂ N ∂ Nη u T4 = V F · C jk : ⊗ ui · F jk · · δu j 2 ∂X ∂X j=1 i=0 (78) ( ) k Nj N jki kk j ∑∑ 1 k jk ∂ N ∂ N T η u + V F · C jk : F jk · ui ⊗ · · δu j . 2 ∂X ∂X j=1 i=0

20

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

Note, that the inner sum operates on the family of particle j from the outer sum. Like before, a rearrangement in index notation merges the two occurring terms: ) ] [ ( Nk ∑ Nj Nk ∑ Nj jki kk j ∑ ∑ ∂ N ∂ N T η u ij j k jk jk jk i T4 = ·C · ·F ·u = δu · V F · δu j · M2 · ui , (79) ∂X ∂X j=1 i=0 j=1 i=0 with the second substantial matrix ) ( jki kk j ∂ Nη T ij jk ∂ Nu k jk ·C · · F jk M2 = V F · ∂X ∂X or in Voigt notation ( ) T ij jT · D jk · Bujki · F jk . M2 = V k F jk · Bkk η

(80)

(81)

References [1] S.A. Silling, Reformulation of elasticity theory for discontinuities and long-range forces, J. Mech. Phys. Solids 48 (1) (2000) 175–209. [2] S.A. Silling, Dynamic fracture modeling with a meshfree peridynamic code, in: Computational Fluid and Solid Mechanics 2003, Elsevier, 2003, pp. 641–644. [3] A. Javili, R. Morasata, E. Oterkus, S. Oterkus, Peridynamics review, Math. Mech. Solids (2018). [4] S.A. Silling, M. Epton, O. Weckner, J. Xu, E. Askari, Peridynamic states and constitutive modeling, J. Elasticity 88 (2) (2007) 151–184. [5] T.L. Warren, S.A. Silling, A. Askari, O. Weckner, M.A. Epton, J. Xu, A non-ordinary state-based peridynamic method to model solid material deformation and fracture, Int. J. Solids Struct. 46 (5) (2009) 1186–1195. [6] M.A. Bessa, J.T. Foster, T. Belytschko, W.K. Liu, A meshfree unification: reproducing kernel peridynamics, Comput. Mech. 53 (6) (2014) 1251–1264. [7] W.K. Liu, S. Jun, S. Li, J. Adee, T. Belytschko, Reproducing kernel particle methods for structural dynamics, Internat. J. Numer. Methods Engrg. 38 (10) (1995) 1655–1679. [8] G.C. Ganzenmüller, S. Hiermaier, M. May, On the similarity of meshless discretizations of peridynamics and smooth-particle hydrodynamics, Comput. Struct. 150 (2015) 71–78. [9] L.B. Lucy, A numerical approach to the testing of the fission hypothesis, Astron. J. 82 (1977) 1013–1024. [10] R.A. Gingold, J.J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (3) (1977) 375–389. [11] J.J. Monaghan, SPH without a tensile instability, J. Comput. Phys. 159 (2) (2000) 290–311. [12] Stewart A. Silling, Introduction to peridynamics, in: Handbook of Peridynamic Modeling, Chapman and Hall/CRC, 2016, pp. 63–98. [13] John T. Foster, Constitutive modeling in peridynamics, in: Handbook of Peridynamic Modeling, Chapman and Hall/CRC, 2016, pp. 181–216. [14] D.J. Littlewood, Simulation of dynamic fracture using peridynamics, finite element modeling, and contact, in: ASME 2010 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers, 2010, pp. 209–217. [15] M.S. Breitenfeld, P.H. Geubelle, O. Weckner, S.A. Silling, Non-ordinary state-based peridynamic analysis of stationary crack problems, Comput. Methods Appl. Mech. Engrg. 272 (2014) 233–250. [16] S.A. Silling, Stability of peridynamic correspondence material models and their particle discretizations, Comput. Methods Appl. Mech. Engrg. 322 (2017) 42–57. [17] C.T. Wu, B. Ren, A stabilized non-ordinary state-based peridynamics for the nonlocal ductile material failure analysis in metal machining process, Comput. Methods Appl. Mech. Engrg. 291 (2015) 197–215. [18] A. Yaghoobi, M.G. Chorzepa, Higher-order approximation to suppress the zero-energy mode in non-ordinary state-based peridynamics, Comput. Struct. 188 (2017) 63–79. [19] X. Gu, E. Madenci, Q. Zhang, Revisit of non-ordinary state-based peridynamics, Eng. Fract. Mech. (2017). [20] E. Madenci, A. Barut, M. Futch, Peridynamic differential operator and its applications, Comput. Methods Appl. Mech. Engrg. 304 (2016) 408–451. [21] S.R. Chowdhury, P. Roy, D. Roy, J.N. Reddy, A modified peridynamics correspondence principle: Removal of zero-energy deformation and other implications, Comput. Methods Appl. Mech. Engrg. 346 (2019) 530–549. [22] E. Madenci, E. Oterkus, Peridynamic Theory and its Applications, Springer, 2016. [23] P. Wriggers, Nonlinear Finite Element Methods, Springer Science & Business Media, 2008. [24] Y. Krongauz, T. Belytschko, Consistent pseudo-derivatives in meshless methods, Comput. Methods Appl. Mech. Engrg. 146 (3–4) (1997) 371–386. [25] P. Lancaster, K. Salkauskas, Curve and Surface Fitting: An Introduction, Academic Press, 1986. [26] M. Arroyo, M. Ortiz, Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods, Internat. J. Numer. Methods Engrg. 65 (13) (2006) 2167–2202. [27] C. Weißenfels, P. Wriggers, Stabilization algorithm for the optimal transportation meshfree approximation scheme, Comput. Methods Appl. Mech. Engrg. 329 (2018) 421–443. [28] J. Korelc, P. Wriggers, Automation of Finite Element Methods, Springer, 2016.

T. Bode, C. Weißenfels and P. Wriggers / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112636

21

[29] M.D. Brothers, J.T. Foster, H.R. Millwater, A comparison of different methods for calculating tangent-stiffness matrices in a massively parallel computational peridynamics code, Comput. Methods Appl. Mech. Engrg. 279 (2014) 247–267. [30] D.J. Littlewood, Roadmap for peridynamic software implementation, SAND Report, Sandia National Laboratories, Albuquerque, NM and Livermore, CA, 2015. [31] C. Kadapa, W.G. Dettmer, D. Peri´c, Subdivision based mixed methods for isogeometric analysis of linear and nonlinear nearly incompressible materials, Comput. Methods Appl. Mech. Engrg. 305 (2016) 241–270. [32] G.C. Ganzenmüller, An hourglass control algorithm for Lagrangian smooth particle hydrodynamics, Comput. Methods Appl. Mech. Engrg. 286 (2015) 87–106.