A semi-Lagrangian constitutive correspondence framework for peridynamics

A semi-Lagrangian constitutive correspondence framework for peridynamics

Journal Pre-proof

A semi-Lagrangian constitutive correspondence framework for peridynamics Masoud Behzadinasab, John T. Foster PII: DOI: Reference:

S0022-5096(19)30951-2 https://doi.org/10.1016/j.jmps.2019.103862 MPS 103862

To appear in:

Journal of the Mechanics and Physics of Solids

Received date: Revised date: Accepted date:

8 October 2019 20 November 2019 28 December 2019

Please cite this article as: Masoud Behzadinasab, John T. Foster, A semi-Lagrangian constitutive correspondence framework for peridynamics, Journal of the Mechanics and Physics of Solids (2020), doi: https://doi.org/10.1016/j.jmps.2019.103862

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

A semi-Lagrangian constitutive correspondence framework for peridynamics Masoud Behzadinasaba,∗, John T. Fostera,b,∗ a Department b Hildebrand

of Aerospace Engineering & Engineering Mechanics, The University of Texas at Austin, United States Department of Petroleum & Geosystems Engineering, The University of Texas at Austin, United States

Abstract Lagrangian formulations, in which material models are established upon a mapping from a deformed domain back to an undeformed reference configuration, have traditionally been employed in most previous developments of the peridynamic theory. In this paper, a semi-Lagrangian constitutive framework for peridynamics is proposed, in which peridynamic material point interactions depend only on their current properties (e.g. position and stress values) in the deformed configuration. A nonlocal version of the velocity gradient is proposed to determine the Cauchy stress rate, using local constitutive theories, as an intermediate quantity in computing peridynamic bond forces. It is shown that a pure state-based peridynamic approach would involve material instabilities. A bond-associated formulation is presented to benefit from both bond-based and state-based views, resulting in an enhanced stability. A bond-associated, correspondence damage modeling is introduced to incorporate classical failure criteria in a convenient way. The proposed framework opens up new pathways in modeling extreme events, where excessive material deformation and damage make the Lagrangian formulation unsuitable. Capabilities of the new theory is showcased in a handful of illustrative scenarios. Keywords: Peridynamics, semi-Lagrangian, correspondence modeling, extreme events, nonlocal

1. Introduction

5

10

15

20

25

Peridynamics is a nonlocal theory of the mechanics of solid deformations. It was originated by Silling (2000) primarily targeting a mathematical consistent representation of the mechanics of continuous and discontinuous media, which was achieved by replacing the partial spatial derivatives in the governing equations with an integral formulation. Given its natural capabilities in modeling material damage without the need for complex numerical treatments at discontinuities, i.e. cracks, the theory has gained broad interest from the computational fracture mechanics community. Numerous constitutive theories have been proposed for peridynamics since its introduction (Silling et al., 2007; Foster et al., 2010; Mitchell, 2011; Tupek and Radovitzky, 2014; Foster and Xu, 2018). Peridynamics has been used to model various challenging problems, e.g. modeling nanostructures (Silling and Bobaru, 2005), fracture in concrete (Gerstle et al., 2007; Yaghoobi and Chorzepa, 2015; Nikravesh and Gerstle, 2018), crack branching phenomenon (Ha and Bobaru, 2010), thermo-mechanical failure analysis in electronic packages (Oterkus et al., 2014), corrosion damage (Jafarzadeh et al., 2018), intra-ply fracture in composites (Mehrmashhadi et al., 2019), and meso-scale modeling of intra-granular fracture (Behzadinasab et al., 2018). Peridynamics has also been coupled with local theories, as well as other meshless methods, to develop efficient tools for failure analysis of engineering structures (Kilic and Madenci, 2010; Seleson et al., 2015; Shojaei et al., 2016; Zaccariotto et al., 2018). There is a special class of peridynamic material models, called constitutive correspondence models, which provide a convenient approach for incorporating classical constitutive models within the peridynamic framework. In peridynamic correspondence models a kinematic variable is utilized as an intermediate quantity in calculating the material point interaction potentials. Originally, Silling et al. (2007) introduced a nonlocal deformation gradient tensor to adapt classical material models through its work conjugate, the first Piola-Kirchhoff stress tensor. More recently, Foster and Xu (2018) presented a finite deformation correspondence theory and used a nonlocal family of Seth-Hill strain measures as the intermediate quantity in computing the peridynamic bond forces. Despite offering convenience, material instabilities have been associated with correspondence material models both within their continuum descriptions and the discretizations (Littlewood, 2011; Tupek and Radovitzky, 2014; Silling, 2017). Zero-energy mode oscillations appearing in meshless peridynamic simulations have been primarily attributed to the integration averaging that allows cancellation of nonuniform parts of deformation and that leads to a drastic degradation in the accuracy of the predicted displacement fields, particularly in regions of steep gradients (Breitenfeld et al., 2014; Tupek and Radovitzky, 2014; Silling, 2017; Chowdhury et al., 2019). Another type of material instability was recently reported to occur even for homogeneous deformations within the finite deformation correspondence framework (Behzadinasab and Foster, 2020). The issue is largely related to a weak dependence of each material point interaction on its own deformation because of the large ∗ Corresponding

authors. Email addresses: [email protected] (Masoud Behzadinasab), [email protected] (John T. Foster)

Preprint submitted to Journal of the Mechanics and Physics of Solids

January 2, 2020

30

35

40

45

50

55

60

influence of the nonlocal integral operator. Several attempts have been recently proposed to improve the stability of correspondence materials. Littlewood (2010); Silling (2017) added penalty terms to the governing equations, which resists deviation from uniform deformations. Chen (2018) and Breitzman and Dayal (2018) reduced the average nodal integration impact and increased the influence of each bond deformation in evolving its own force state. An arbitrary split of a neighborhood of peridynamic bonds into multiple sub-regions (Chowdhury et al., 2019) enhanced the stability in some problems, but is unsettling as a continuum theory. Most developments of the peridynamic theory have taken place using a Lagrangian formulation, in which the material point interactions depend on their current positions with respect to an undeformed reference configuration. The Lagrangian framework can become unreliable once very large deformations occur. Bergel and Li (2016) proposed an updated-Lagrangian peridynamic model that forms peridynamic neighborhoods based on the deformed configuration and describes the relative degree of interaction between material points (aka influence function) using the current positions. However, their theory still involves a mapping between the deformed and undeformed domains, which can become severely distorted under extremely large deformations. Silling et al. (2017) were the first to present an Eulerian peridynamic model that defines bond forces based on only the current configuration. They showed thermodynamic consistency of the formulation and demonstrated the applicability of the theory to simulating shock waves and fluid motion. While the Eulerian model has certain applications, it is difficult to simulate historydependent problems, e.g. cyclic loading and metal plasticity, within the framework. Regardless, the approach of Silling et al. (2017) does not have the convenience of correspondence material models, and a complete recast of a classical constitutive model is required to recover analogue physical behavior. In addition, surface effects, a well-known issue in peridynamics (Mitchell et al., 2015), can be a point of concern within the Eulerian theory. In this work, a semi-Lagrangian constitutive correspondence framework for peridynamics is developed. Similar to the Eulerian model, only the positions of material points in the current deformed configuration determine their relative degree of interaction. The stress rate measures only depend on the current configuration of material points. A nonlocal analog of the spatial velocity gradient tensor is proposed to evolve the stress measure and the corresponding bond forces. However, internal state variable are stored at Lagrangian material points, thus enabling history dependence, in contrast to a full Eulerian formulation. The presented model does not suffer from any surface effects and can be used to incorporate classical constitutive equations in a straightforward manner. The remainder of the paper is presented as follows: Section 2 provides an overview of the peridynamic material modeling. In Section 3 a semi-Lagrangian peridynamic framework is introduced, and energy balance within its context is studied. Section 4 presents a correspondence material model for semi-Lagrangian peridynamics and identifies a material instability problem if a completely state-based approach is taken. In Section 5 the initial formulation is modified, and a stable, bond-associated counterpart is introduced. In Section 6 a correspondence damage model is given for the bond-associated theory. Capabilities of the semi-Lagrangian framework is showcased through numerical examples in Section 7. Section 8 provides the concluding remarks. 2. Peridynamics background The classical (local) continuum theory of the momentum conservation for solid mechanics in Lagrangian form is stated as ¨ ρ0 u(X, t) = ∇ · P(X, t) + ρ0 b(X, t),

65

(1)

in which ρ0 is the undeformed mass density, X is a material point in the undeformed reference configuration, u is the displacement field, and b is an external body force. The first term on the right-hand side of Eq. (1) is the divergence of the first Piola-Kirchhoff stress P that is governed by material deformation. Therefore, classical mechanics involves the use of spatial derivatives of displacements, which are undefined at singularity points. Consequently, any discontinuity in the displacement field (e.g. cracks) is inconsistent with the standard classical formulation. While special numerical techniques, e.g. XFEM (Moës et al., 1999), have been proposed to simulate cracks, they involve additional numerical complexities and supplemental governing equations for determining the path of crack growth. On the other hand, the nonlocal theory of peridynamics has been presented (Silling, 2000; Silling et al., 2007) in a way that embraces any material discontinuity. The Lagrangian peridynamic model for linear momentum balance is Z   ¨ ρ0 u(X, t) = T(X, t)hξi − T(X + ξ, t)h−ξi dξ + ρ0 b(X, t), H(X)

which replaces the partial differential equation form of the classical equation Eq. (1) with an integral functional that is used to describe the constitutive interactions of material points through a vector field T (with units of force per square volume). H(X) is the family set of X defined in the reference configuration by H(X) = {X + ξ | X + ξ ∈ B0 , 0 < |ξ| ≤ δ} ,

70

in which the vector ξ is named a peridynamic bond. The characteristic length scale, typically the radius of a spherical neighborhood is defined with δ, called the horizon. The material point X is excluded from its own family. An example of a peridynamic 2

body is shown in Fig. 1.

φ B0

Bt ξ x(X, t)hξi

X

x(X, t)

H(X)

Ht (X) Figure 1: Example of a peridynamic body mapped by a deformation state from its reference configuration to a deformed domain.

The force vector state T is the Fréchet derivative of a strain-energy functional ψ (Foster and Xu, 2018, Section 2), i.e., Thξi = ∇ψ(xhξi), where ψ is assumed to be a function of the deformation vector state x, i.e., x = x(X, t)hξi = u(X + ξ, t) − u(X, t) + ξ,

= x(X + ξ, t) − x(X, t),

75

where x(X, t) is the position of material point X at time t. Note that in the peridynamic literature the deformation state is usually denoted by Y (where y is used to denote material positions in the current domain). A peridynamic vector state is a nonlinear operator that maps from H to R3 (a scalar state is a mapping to R). For example, x : V → V ∈ R3 transforms all vectors ξ originating from X in the reference configuration into their deformed images x in the current configuration. As the displacement field u can be discontinuous in general, the deformation state does not need to be a linear operator and, in general, exhibits a relaxation in kinematics compared with classical mechanics. 2.1. Material modeling

80

85

90

95

The constitutive behavior of the material is governed by the specific choice of the stored energy density ψ. Peridynamic material models are generally categorized via different classes: bond-based vs. state-based, and ordinary vs. non-ordinary. The original peridynamic theory (Silling, 2000) is a bond-based formulation, which determines a bond force solely based on the deformation history of each bond independent of all other conditions. To address the shortcomings associated with the bondbased approach, e.g. limitation to an effective Poisson ratio of 1/4 in three dimensions and inconsistency with metal plasticity, Silling et al. (2007) proposed a generalized form of peridynamics, called the state-based theory, which allows the behavior of a material point to depend on the collective deformation history of all bonds within H. In ordinary material models, force vector states are parallel to the deformation vector states. Non-ordinary materials have force vector states that are not parallel to the deformation vector states. Ordinary materials are automatically equipped with angular momentum conservation, while additional restrictions are required to maintain this conservation property in non-ordinary materials (Silling et al., 2007). 2.1.1. Correspondence materials Silling et al. (2007) proposed the notion of constitutive correspondence which involves using a nonlocal analog of a kinematic variable (e.g. the deformation gradient) as a connection between the local (classical) and nonlocal peridynamic frameworks. Therefore, correspondence materials do not require a complete recast of classical material models for replicating their constitutive behaviors and offer a notable convenience for incorporating local models within the peridynamic setting. A nonlocal version of the deformation gradient tensor is defined as (Silling et al., 2007) ! Z ¯F = ωhξi xhξi ⊗ ξ dξ K−1 , (2) H

where K is a symmetric second-order shape tensor, given by Z K= ωhξi ξ ⊗ ξ dξ. H

3

The relative degree of interactions between the material points are described by the scalar state ωhξi, called the influence state. It was shown that under a continuous, homogeneous deformation the nonlocal deformation gradient recovers exactly the local deformation gradient (Silling et al., 2007, Section 12). A Lagrangian local theory defines a strain-energy functional based on the work conjugate of the deformation gradient tensor and the first Piola-Kirchhoff stress tensor. Therefore, the nonlocal kinematic ¯ which results in the following force vector state quantity can be used as a metric to incorporate local theories using P = P(F), (Silling et al., 2007, Section 18) −1 T k = ωhξi Pkp K pq ξq . Additionally, a finite deformation constitutive correspondence has been presented (Foster and Xu, 2018) using a different kinematic variable. A family of nonlocal right Cauchy-Green deformation tensors is defined as ! Z  xhξi 2m ξ ⊗ ξ   ¯ ωhξi C(m) =  dξ : L−1 , m , 0, |ξ| |ξ|2 H

where L is a symmetric fourth-order shape tensor, given by Z ξ⊗ξ⊗ξ⊗ξ L= ωhξi dξ. |ξ|4 H

A set of nonlocal Lagrangian Seth–Hill strains is then introduced as  1 ¯ C(m) − I , m , 0, E¯ (m) = 2m and !  Z xhξi ξ ⊗ ξ   ¯  dξ : L−1 , E(0) =  ωhξi ln |ξ| |ξ|2 H

(3)

which is shown to retrieve exactly their local counterparts under homogeneous setting (Foster and Xu, 2018, Section 3). Employing a local hyperelastic constitutive law within the peridynamic model, the force vector state is resulted as ! ξ p ξq −1 |x| 2m xk T k = ωhξi S (m)i j 2 L pqi j , (4) |ξ| |ξ| |x|2 ¯ (m) ). where the generalized Kirchhoff stress S(m) is the energy conjugate of E¯ (m) , i.e., S(m) = S(m) (E 2.2. Fracture modeling

100

105

110

In peridynamics, material failure is typically achieved through irreversible bond breakage. The critical stretch criterion (Silling and Askari, 2005) is a widely-used bond-based damage model in the peridynamic literature. In this model, a peridynamic bond is broken once its stretch, i.e. the ratio of the deformed length to the initial length, exceeds a prescribed value, called the critical stretch which can be related to the critical energy release rate Gc . Once a bond breaks, the bond force is redistributed to the other (unbroken) bonds, leading to damage coalescence and progressive crack propagation. The critical stretch failure model only considers the volumetric part of deformation and does not engage the deviatoric part. Despite certain applications to brittle solids, the damage theory is not suitable for modeling fracture in metallic alloys where elasto-plastic behavior can play a large role. To address this shortcoming, Foster et al. (2011) presented an energy-based failure criterion where a bond breaks after it reaches a certain critical work performed on it through deformation. The critical energy density level for failure can be determined form the critical energy release rate. A state-based damage correspondence framework was proposed by Tupek et al. (2013) to integrate local failure theories conveniently into the peridynamic formulation. In this model, local conditions, i.e. plasticity, strain rate, stress triaxiality, temperature, etc., evolve the state of a material point damage through a (classical) continuum damage model. After a material point is damaged, all of its interactions with the neighboring points are degraded (through decaying the influence states), and the force carried by its bonds are transferred to the other bonds. 3. Semi-Lagrangian peridynamics

115

With the motivation of modeling extreme events, i.e. significantly large deformations, a semi-Lagrangian peridynamic formulation is introduced in this section. In this model material points are still tracked in motion (the Lagrangian view), rather than looking at fixed points in the space (the Eulerian view). However, an updated kernel defines the degree of interactions between material points. Fig. 2 compares a semi-Lagrangian kernel with a Lagrangian one. A Lagrangian kernel may become completely distorted in very large deformations, and the Lagrangian kinematic variables defined in Eqs. (2) and (3) can become non-physical, i.e., have non-positive determinants.

4

φ Lagrangian scheme

X

x(X, t)

H(X, t)

H0 (X)

φ Semi-Lagrangian scheme

X

x(X, t)

H(X, t)

H0 (X)

Figure 2: Lagrangian versus semi-Lagrangian discrete schemes. In semi-Lagrangian, the kernel preserves its shape regardless of deformation.

120

3.1. Governing equations Let x describe position of a body B in its current deformed configuration. Suppose ρ = ρ(x) is the current density field, and b = b(x) is an applied external body force field. The internal force density is described with f = f(x). The balance of linear momentum on a subregion P ⊂ B reads Z Z Z d ρ˙x dx = f dx + ρb dx. dt P P P Invoking the mass balance principle,

Z

P

where

ρ

D˙x dx = Dt

Z

f dx +

P

Z

ρb dx,

(5)

P

D() is the material derivative operator. Rearranging Eq. (5), Dt Z   D˙x ρ − f − ρb dx = 0. Dt P

As the subregion P is arbitrary, the local balance of linear momentum is resulted as ρ

D˙x = f + ρb. Dt

(6)

Remark 1. Due to the material derivative operator, the Lagrangian approach is the preferred reference frame to employ Eq. (6) in modeling solid mechanics. By expanding the material derivative, the balance equation is ρ

∂˙x = f − ρ˙x · L + ρb, ∂t

(7)

∂˙x where L = is the spatial velocity gradient tensor. A nonlocal version of this kinematic variable is introduced later (see ∂x Sections 4 and 5). An Eulerian approach is the preferred reference frame for representing the local linear momentum balance as in Eq. (7) with particular utility in fluid mechanics. Similar to the Lagrangian state-based peridynamic theory, we define the internal state of a material point as the summation

5

of its collective interactions with other material points in the body Z h i f(x) = t(x)hηi − t(x + η)h−ηi dη,

(8)

H(x)

where H(x) contains all the neighboring points with which x directly interacts in the current domain, i.e.,  H(x) = x + η | x + η ∈ B, 0 < |η| ≤ δ . 125

η refers to a peridynamic bond in the deformed frame. t is still named the force vector state, but in the current configuration. Note that t can be seen as similar to the Cauchy stress (or true stress) and T to the first Piola-Kirchhoff stress (or engineering stress), considering the frames in which they operate. There is a tensor-valued functional Λ in local hypoelastic theories such that ˙ = Λ(σ, L), σ

˙ is the material rate of Cauchy stress tensor. Similarly, the force state can be generally expressed as where σ ˙t = Λ(t, L, x˙ ), where x is the deformation vector state in the context of the semi-Lagrangian model, i.e., x = x(x)hηi = η. A material behavior can be path dependent, in general, governed by a flow equation. If there exists a strain-energy functional describing the response of matter, the material is hyperelastic, i.e., t = Γ(x). 3.2. Energy balance An expression of energy balance can be obtained directly from the balance of linear momentum, i.e., Eq. (6), by multiplying by x˙ and integrating over an arbitrary subregion P Z Z Z d ρ˙x · x˙ dx = F · x˙ dx + ρb · x˙ dx. (9) dt P P P

Using Eqs. (8) and (9),

d dt

Z

P

ρ˙x · x˙ dx =

Z Z P

H(x)

Z h i t(x)hηi − t(x + η)h−ηi · x˙ dη dx + ρb · x˙ dx.

For brevity, use x0 = x + η, t = t(x)hηi, and t0 = t(x + η)h−ηi. Then, Z Z Z Z Z Z     0 0 0 0 t · x˙ − t · x˙ dx dx − t − t · x˙ dη dx = It is trivial that

P

Using Eqs. (11) and (12), Z Z P

P

H(x)

H(x)

H(x)

Z Z P



H(x)∩P

P



P

H(x)

 t · x˙ 0 − t0 · x˙ dx0 dx = 0.

Z Z  t − t0 · x˙ dη dx =

H(x)\P

(10)

P


t · x˙ 0 − x˙ dx0 dx.

(11)

(12)

Z Z   t · x˙ 0 − t0 · x˙ dx0 dx − P

H(x)


t · x˙ 0 − x˙ dx0 dx. | {z } η˙

(13)

Substituting Eq. (13) in Eq. (10), Z Z Z Z Z Z   d 0 0 0 ρ˙x · x˙ dx + t · η˙ dη dx = ρb · x˙ dx. t · x˙ − t · x˙ dx dx + dt P P H(x) P H(x)\P P

The first term on the left-hand side is the rate of change of kinetic energy of P. The first term on the right side of the equation is the contributions of the region outside of P, thus a form of supplied power, similar to the external body force role. Therefore, Z Z ˙ ˙ sup (P). K(P) + t(x)hηi · η˙ dη dx = W (14) P

H(x)

Comparing Eq. (14) with the first law of thermodynamics, which theorize existence of an internal energy U such that ˙ ˙ ˙ sup (P), K(P) + U(P) =W where heat transfer is ignored, results in ˙ U(P) =

Z Z P

H(x)

t(x)hηi · η˙ dη dx.

6

(15)

Setting P = B, the global form of Eq. (15) is obtained as Z Z ˙ U(B) = B

H(x)

t(x)hηi · η˙ dη dx.

Note that if t(x)hηi = 0 for |η| > δ, i.e., if x0 is not a neighbor of x, Eq. (16) is equivalent to Z Z
˙ U(B) = t(x)hx0 − xi · x˙ 0 − x˙ dx0 dx. B

(16)

(17)

B

For more discussions on energy balance, including heat generation, refer to (Silling and Lehoucq, 2010, Section 2.5).

4. Correspondence material modeling 130

A constitutive correspondence model for the semi-Lagrangian formulation is developed in this section. A nonlocal velocity gradient operator is used as the intermediate quantity in incorporating local material models to evaluate the force state. 4.1. Peridynamic spatial derivative operator Bergel and Li (2016) proposed a peridynamic spatial derivative operator as ! Z D{g} = ωhηi ζ(g)hηi ⊗ η dη M−1 ,

where g = g(x) is a given function, ζ(g) is defined as

(18)

H

ζ(g)hηi = ζ(g)[x]hηi = g(x + η) − g(x), ω is an influence function, which is set to 0 if the relative distance between material points is greater than the horizon, i.e., ωhηi = 0, M is called the spatial shape tensor, given as M=

Z

H

135

140

if |η| > δ.

ωhηi η ⊗ η dη.

(19)

The use of a shape tensor, similar to the Lagrangian correspondence theory, eliminates or minimizes the surface effects. The nonlocal operator D recovers its local counterpart, i.e. the local spatial derivative operator, under homogeneous deformations or as the horizon approaches zero (Bergel and Li, 2016). Remark 2. In practice, as the horizon is finite, D approximates the zeroth-order local spatial derivative operator (zeroth-order ∂g in the sense that must be uniform throughout the field – g can vary linearly with x). For obtaining higher-order nonlocal ∂x analogs of the local quantity, a similar approach to the implicit gradient model (Chen et al., 2004) from reproducing kernel particle methods (RKPM) (Liu et al., 1995) or the peridynamic (Lagrangian) differential operator (Madenci et al., 2016) can be followed, with its adaptation to the current configuration. That is not our focus in this work, and we suffice to the zeroth-order. Utilizing a higher-order nonlocal derivative operator in this framework, however, must be straightforward. 4.2. State-based velocity gradient Using the velocity field as g in the peridynamic spatial derivative operator, i.e. Eq. (18), a nonlocal velocity gradient is obtained: ! Z ¯ = L ωhηi η˙ ⊗ η dη M−1 , (20) H

¯ is called the state-based velocity gradient since in order to compute velocity gradient for a material point, an average over all L its associated bonds takes place.

Remark 3. The velocity gradient operator can be used to update the volumes in the current configuration. A deformed volume element dx is related to its reference state dX by dx(t) = J(X, t) dX, where J is the Jacobian determinant of the deformation gradient in the classical theory, in general it describes the change in differential volumes between the two configurations, and can be evolved through the following integration Z t ˙ dτ, ˙ = J(τ) L¯ kk (τ). J(t) = J(τ) J(τ) 0

7

Remark 4. In applications of the Lagrangian correspondence theory, where a nonlocal deformation gradient F¯ is employed, the velocity gradient is computed using L¯ = F˙¯ F¯ −1 145

150

to calculate stress increments. For example, the Flanagan-Taylor algorithm (Flanagan and Taylor, 1987) has been employed in correspondence material simulations (Warren et al., 2009; Foster et al., 2010; Amani et al., 2016) to compute the polar decomposition efficiently to "un-rotate" the Cauchy stress rate for a co-rotational, i.e. objective constitutive modeling theory. Having the ability to directly evaluate the velocity gradient tensor, as in Eq. (20), can improve accuracy. 4.3. Force vector state evaluation In order to incorporate local material models, let us assume there exists a strain-energy functional that describes the behavior of a peridynamic body, i.e., the material is hyperelastic. The force state can then be evaluated using Eq. (16). The classical theory gives the power conjugate of a deforming hyperelastic body B as Z ˙ U(B) = σ : L dx. (21) B

Using the proposed peridynamic spatial velocity gradient tensor in Eq. (20), let’s introduce a peridynamic hyperelastic material such that Z ˙ U(B) = σ : L¯ dx, ZB = σi j Li j dx, B ! Z Z ωhηi η˙ i η p dη M −1 = σi j (22) p j (x) dx. B

H(x)

Comparing Eq. (16) and Eq. (22), and using the symmetry in σ,

ti = ωhηi η p M −1 p j σ ji , or t(x)hηi = ωhηi η M−1 (x) σ(x).

(23)

Returning to the classical theory, the Cauchy stress tensor is governed through ∂˙u , (24) ¯ ∂L where u˙ is a strain power density function (in the deformed domain). The peridynamic constitutive correspondence model is completed using Eqs. (20), (23) and (24) together. It is realized from Eq. (23) that the force state is not (necessarily) in the deformed bond direction. Therefore, the material model is non-ordinary, and the balance of angular momentum is not automatically satisfied and needs to be analyzed. Let i jk denote the Levi-Civita symbol. Then the net moment of the force state t on a material point is calculated using Eq. (23) and the symmetry of σ: ! Z Z η × thηi dη = i jk ηi σ jl ωhηi η p M −1 pl dη, H H i Z = i jk σ jl M −1 ωhηi ηi η p dη, pl σ=

H

= i jk σ jl M −1 pl Mip ,

= i jk σ jl δil , = i jk σ ji , = 0, thus, the balance of angular momentum holds.

155

Remark 5. Matter interpenetration can happen in state-based correspondence materials since the effect of an individual bond deformation is damped through the (averaging) integral operation (the collective deformation of neighborhood determines the stress) (Tupek and Radovitzky, 2014). In the context of the semi-Lagrangian framework, collapse of a sub-horizon into a single point can take place while the nonlocal velocity gradient and, subsequently, the Cauchy stress measure remain finite. The interested reader can refer to (Tupek and Radovitzky, 2014, Section 3) for more discussions on kinematic limitations of (Lagrangian) correspondence materials. The issue can be naturally handled in the semi-Lagrangian numerical simulations by carefully choosing an influence function that results in infinite bond forces when two peridynamic nodes are passing through each other in a meshless particle setting, i.e., 8

η → 0 ⇒ |thηi| → ∞.

(25)

Suppose ωhηi → |η|−m as η → 0. Also, let σ and M have finite values. Referring back to Eq. (23), η → 0 ⇒ |thηi| → |η|−m+1 .

It is immediately clear that choosing m > 1 satisfies Eq. (25), therefore automatically enforces the non-interpenetration of matter. Note that having m ≥ 3 would lead to an ill-conditioned shape tensor in continuum scale. Our favorite choice for the influence function is ! |η| 1 ωhηi = 2 ω ˜ , (26) δ |η| 160

with a condition on ω ˜ such that ω ˜ (β) → c > 0 as β → 0. The definition in Eq. (26) results in a shape tensor with units of volume. In addition, using Eq. (20), it is seen that this choice of influence function would further secure ||L|| → ∞ as η → 0 for some bonds in the neighborhood. 4.4. Material stability

165

Silling (2017) introduced a technique for investigating the material stability of a peridynamic constitutive model using an energy minimization approach, which is equivalent to the Lyapunov stability test (Behzadinasab and Foster, 2020). A similar tool for the presented semi-Lagrangian formulation is developed here. The total potential energy of a bounded elastic peridynamic body B is defined as Π(B) = U(B) − Wb (B),

(27)

where U(B) is the total strain energy and Wb (B) is the total work of the external body force. A given deformation, which satisfies the essential boundary conditions, is an equilibrium if Π(B) is stationary, i.e., δΠ(B) = 0

(28)

for all admissible variations δu. Using Eq. (17) directly Z Z
δU(B) = t(x)hx0 − xi · δu(x0 ) − δu(x) dx0 dx, ZB ZB   = t(x0 )hx − x0 i − t(x)hx0 − xi · δu(x) dx0 dx. B

Using Eqs. (27) and (29) and

δWb (B) = δΠ(B) = −

Z

B

(29)

B

Z

B

ρ(x) b(x) · δu(x) dx,

! Z   0 0 0 0 t(x)hx − xi − t(x )hx − x i dx + ρ(x) b(x) · δu(x) dx.

(30)

B

As δu is arbitrary, for Eq. (28) to hold Z   t(x)hx0 − xi − t(x0 )hx − x0 i dx0 + ρ(x) b(x) = 0, B

∀x ∈ B,

(31)

which is the peridynamic equation of motion, i.e., Eq. (6), in the static form. An equilibrium point is stable if the behavior of the potential energy functional near its vicinity is characterized as convex, i.e., δ2 Π(B) > 0, δΠ(B) = 0,

meaning that any departure from the equilibrium state increases the overall potential energy. Using Eq. (30) and noting Eq. (31), Z Z   2 δ Π(B) = − δt(x)hx0 − xi − δt(x0 )hx − x0 i dx0 · δu(x) dx, B B Z Z
= δt(x)hx0 − xi · δu(x0 ) − δu(x) dx0 dx > 0, B

or

B

δ2 Π(B) =

Z Z B

B

δt(x)hx0 − xi · δη dη dx > 0,

where δη = δu(x0 ) − δu(x). A sufficient condition for holding Eq. (32) is Z δt(x)hx0 − xi · δη dη > 0, ∀x ∈ B and all admissible δη. B

9

(32)

(33)

For the proposed state-based correspondence model, using Eq. (23), δt(x)hηi = δωhηi η M−1 (x) σ(x) + ωhηi δη M−1 (x) σ(x) + ωhηi η δM−1 (x) σ(x) + ωhηi η M−1 (x) δσ(x), thus from Eqs. (33) and (34), if # Z "  δωhηi η M−1 (x) σ(x) + ωhηi δη M−1 (x) σ(x) + ωhηi η δM−1 (x) σ(x) + ωhηi η M−1 (x) δσ(x) · δη dη > 0

(34)

(35)

B

holds for all admissible deformations, the material model is stable unconditionally – note the numerical stability is a different matter, where particle discretization can play a large role. There exist material instabilities for any deformation condition where Eq. (35) is violated. For general deformations, analytical stability analysis may not be easily tractable. Thus, let us focus on homogeneous deformations. Limiting the perturbation to a single material point x0 , i.e., δx0 = ε and δx = 0, ∀x , x0 . Also, let us pick x0 in the bulk of material (not near the boundaries) so that it has a full neighborhood in its surrounding (the symmetry is used in the following to obtain closed-form continuum results). For the bonds associated with x0 , δηhx0 − x0 i = δx0 − δx0 = −ε.
Considering a spherical influence function ωhηi = ω s |η| ,

η · δη dω s |η| δωhηi = δ|η| = ω0s |η| . d|η| |η| Using Eqs. (19), (36) and (37),

Z 

(36)

(37)



η · δη η ⊗ η + ω s |η| δη ⊗ η + η ⊗ δη dη, |η| ZB  

η·ε η ⊗ η + ω s |η| − ε ⊗ η − η ⊗ ε dη. = − ω0s |η| |η| B

δM(x0 ) =

ω0s |η|

Because the integrand is an odd function, using the symmetry of the area of integration, δM(x0 ) = 0.

(38)

˙ Similarly, using Eqs. (20) and (36), and noting that for a homogeneous deformation η(t) = α(t) η with α being a scalar function, Z  

η · δη ¯ 0) = α(t) η ⊗ η + ω s |η| α(t) δη ⊗ η + α(t) η ⊗ δη dη M−1 (x0 ), δL(x ω0s |η| |η| B Z  

η·ε η ⊗ η + ω s |η| − ε ⊗ η − η ⊗ ε dη M−1 (x0 ), = α(t) − ω0s |η| |η| B = 0,

(39)

−1

where δM (x0 ) = 0 from Eq. (38) is implicitly used. A classical theory governs the evolution of the stress state through ˙ = f(σ, L). Therefore, using Eq. (39), σ δσ(x0 ) = 0. (40) Plugging Eqs. (36), (38) and (40) in Eq. (35), the stability criterion for homogeneous conditions reduces to Z 
η·ε
  ω0s |η| η + ω s |η| ε dη M−1 σ ε > 0, |η| B

(41)

where the uniform nature of the problem is implicit, i.e., M(x) = M and σ(x) = σ. 4.4.1. Uniform hydrostatic deformations For pure hydrostatic loading conditions σ = σ I,

(42)

where σ is the hydrostatic pressure, and I is the identity tensor. Moreover, regardless of the loading condition, the shape tensor M (in the bulk of the material) reduces to M = MI (43) for spherical influence functions, in which M is a positive constant. Let us limit the perturbation ε to a single coordinate direction, say the x-coordinate, i.e., ε = ε ˆi. (44) ˆ Substituting Eqs. (42)–(44) in Eq. (41), and using η = η1 ˆi + η2 ˆj + η3 k, Z  
 
η1 ε  ˆ η1 i + η2 ˆj + η3 kˆ + ω s |η| ε ˆi dη M −1 I σ I ε ˆi > 0. ω0s |η| |η| B 10

Rearranging terms and noting that integral of odd functions is zero for symmetric areas of integration, Z 
η2
  ω0s |η| 1 + ω s |η| dη σ > 0, M −1 ε2 |η| B and, as M −1 ε2 > 0,

Ω σ > 0, where Ω is a constant, given by

Z 

η2 Ω= ω0s |η| 1 + ω s |η| dη, |η| B

which depends on how ω s is defined. If Ω is positive, e.g. if  
 c2 ω s |η| =   0

170

175

180

185

(45)

if 0 < |η| < δ, otherwise,

Eq. (45) is violated when σ < 0, i.e. when loading is pure compression. On the other hand, for negative Ω, e.g. if  
 c2 /|η|3 if 0 < |η| < δ, ω s |η| =   0 otherwise,

pure dilatation loading, i.e. σ > 0, causes instability. If Ω = 0 the model is neutrally stable for any hydrostatic loading condition, which is problematic in practice (for a neutrally stable correspondence material, unphysical oscillations in numerical simulations emerge (Behzadinasab and Foster, 2020)). Remark 6. As shown, depending on the choice of influence function, pure dilatation or compression loading causes instability in the continuum level; therefore, the semi-Lagrangian, state-based correspondence model suffers from a material instability (rather than a numerical issue). Meshless smoothed particle hydrodynamics (SPH) (Gingold and Monaghan, 1977) is another semi-Lagrangian numerical technique that is known to exhibit tensile instability (Swegle et al., 1995). In the conventional SPH the instability manifests itself as a clumping of nodes and formation of artificial voids. A similar behavior is observed from the proposed peridynamic model in simulating loading scenarios where the stability condition, i.e. Eq. (45), is violated (see Section 7.1). Remark 7. Unstable behavior is not limited to pure hydrostatic loading scenarios. For example, the state-based model exhibits spurious modes of deformation in simulating uniaxial tension in a 3D semi-Lagrangian bar (see Section 7.1). For more discussions on the implications of the existing instability refer to the stability analysis of Behzadinasab and Foster (2020) on another peridynamic model with a similar issue. Remark 8. Other types of instability have been identified within the state-based peridynamic models; the state-based, semiLagrangian theory is no exception. For example, the zero-energy mode oscillations associated to non-uniform deformations and particle discretization are discussed in (Breitenfeld et al., 2014; Silling, 2017; Chowdhury et al., 2019). This type of instability is largely attributed to the averaging integral operator that allows the cancellation of nonuniform parts of deformation. A modified approach is proposed next, aiming to improve the model stability by increasing the influence of each bond’s deformation on its own force vector state and reduce the influence of averaging.

190

5. A bond-associated, semi-Lagrangian correspondence model Material instability, i.e. zero-energy modes, is largely present in the state-based correspondence theories of peridynamics (Breitenfeld et al., 2014; Silling, 2017; Chowdhury et al., 2019; Behzadinasab and Foster, 2020). The problem is caused by a high bias imposed through the use of a non-local operator that effectively averages deformations within a family of bonds. The integral operator is not sensitive enough to large variances in the deformation field and leads to a non-unique mapping between the approximate deformation gradient and force vector state (Tupek and Radovitzky, 2014; Silling, 2017). In other words, statebased correspondence materials have a too weak of a dependence of each bond’s force density on its individual deformation. To enhance the stability, we combine the proposed state-based approach with a bond-associated view, inspired by the work of Breitzman and Dayal (2018). A bond-level velocity gradient is then defined by ! ¯ ¯ + η) ¯ ¯ + η) L(x) + L(x L(x) + L(x η L(x)hηi = + η˙ − · η ⊗ 2, (46) 2 2 |η|

¯ where L(x) is the average velocity gradient at x, defined in Eq. (20). Note the symmetry in this formulation, i.e., L(x)hηi = L(x + η)h−ηi, which is an improvement over what Breitzman and Dayal (2018) proposed (cf. Eq. (3.2) in (Breitzman and Dayal, 2018)) when deriving a similar formula, but for a nonlocal deformation gradient tensor. This symmetry implies that all related 11

195

bond properties, e.g. bond stress, would evolve in a symmetric form. This feature will be certainly useful when the damage model is introduced such that when a bond hx0 − xi breaks, its counterpart hx − x0 i also breaks (see Section 6). Remark 9. A pure bond-based velocity gradient can be proposed as LBB hηi =

200

η˙ ⊗ η , |η|2

(47)

whose value is determined only using the individual bond deformation. The bond-based model can lead to many limitations similar to the Lagrangian bond-based peridynamics (e.g. limitation on the effective Poisson ratio). On the other hand, a pure state-based version is ¯ ¯ + η) L(x) + L(x LSB (x)hηi = , (48) 2 where the collective deformation of the neighborhood governs the property of the associated bonds, resulting in damping of the deformation inhomogeneities through the (averaging) integral operation, which can lead to instability issues (see Section 4.4). The bond-level velocity gradient given in Eq. (46) is a combination of Eqs. (47) and (48), in order to benefit from the advantages of both approaches by eliminating the bond-based limitations through the state-based idea, while improving stability using the bond-based approach. In the proposed model, individual bond deformations play a larger role compared to the statebased model in determining their own kinematic variable (and therefore force state), while still using the information from the neighborhood. The term in the parentheses in Eq. (46) describes the deviation of a bond deformation from the collective deformation of neighborhood. If the deformation is uniform, the mentioned term vanishes. Therefore, the bond-associated theory also corresponds to the local theory under homogeneous deformations. Let’s define a peridynamic material with the following strain power density functional: Z u˙ (x) = α(x)hηi u˙ hηi dη,

(49)

B

where u˙ hηi is the bond-level strain power density functional and αhηi is a normalized influence state (with units of volume−1 ), Z ωhηi α(x)hηi = , ω0 (x) = ωhηi dη. ω0 (x) B

A classical (local) theory postulates that the strain power density is the power conjugate of the velocity gradient and the Cauchy stress tensors in the deformed domain, i.e., ˙ u˙ hηi = ψ(Lhηi) = σhηi : Lhηi, (50) 205

where σhηi is the bond-level Cauchy (true) stress, a quantity depending on the history of Lhηi (see Eq. (24)). Combining Eqs. (20), (46), (49) and (50), the total strain power of the peridynamic body B is obtained as Z ˙ U(B) = u˙ (x) dx, x∈B

=

Z Z

x∈B x0 ∈B

=

Z Z

x∈B x0 ∈B

=

Z Z

x∈B x0 ∈B

=

Z Z

x∈B x0 ∈B

α(x)hx0 − xi u˙ (x)hx0 − xi dx0 dx, α(x)hx0 − xi σi j hx0 − xi Li j hx0 − xi dx0 dx,  ! x0 − x   L¯ i j (x) + L¯ i j (x0 ) j  L¯ ip (x) + L¯ ip (x0 ) 0 ωhx0 − xi j  0 0 0 σi j hx − xi  + x˙i − x˙i − (x p − x p ) 0  dx dx, ω0 (x) 2 2 |x − x|2 ( x˙i0 − x˙i )(x0j − x j ) ωhx0 − xi σi j hx0 − xi dx0 dx ω0 (x) |x0 − x|2

  ! (x0p − x p )(x0j − x j )   −1  dx0 dx + ωhˇx − xi ( x˙ˇi − x˙i )( xˇl − xl ) dˇx Mlp (x) δ p j − |x0 − x|2 x∈B x0 ∈B xˇ ∈B   ! Z Z Z (x0p − x p )(x0j − x j )   ωhx0 − xi 1 −1 0   dx0 dx. + σi j hx0 − xi ωhˇx − x0 i ( x˙ˇi − x˙0 i )( xˇl − xl0 ) dˇx Mlp (x ) δ p j − ω0 (x) 2 |x0 − x|2 Z Z

x∈B x0 ∈B

ωhx0 − xi 1 σi j hx0 − xi ω0 (x) 2

Z

xˇ ∈B

12

Exchanging the limits in the last two terms, making use of the symmetry in bond properties, and rearranging to obtain Z Z ( x˙i0 − x˙i )(x0j − x j ) ωhx0 − xi ˙ U(B) = dx0 dx σi j hx0 − xi ω0 (x) |x0 − x|2 x∈B x0 ∈B ! ! Z Z Z ( xˇ p − x p )( xˇ j − x j ) ωhˇx − xi 1 0 0 0 0 −1 dˇx dx + σi j hˇx − xi ωhx − xi ( x˙i − x˙i )(xl − xl ) dx Mlp (x) δ p j − ω0 (x) 2 |ˇx − x|2 0 x∈B xˇ ∈B x ∈B ! ! Z Z Z ( xˇ p − x p )( xˇ j − x j ) ωhˇx − xi 1 0 0 0 0 −1 + σi j hˇx − xi ωhx − xi ( x˙i − x˙i )(xl − xl ) dx Mlp (x) δ p j − dˇx dx, ω0 (ˇx) 2 |ˇx − x|2 x∈B xˇ ∈B

Z Z

x0 ∈B

( x˙i0 − x˙i )(x0j − x j )

0

ωhx − xi dx0 dx σi j hx0 − xi ω0 (x) |x0 − x|2 x∈B x0 ∈B ! ! # Z Z " Z ( xˇ p − x p )( xˇ j − x j ) 1 1 1 −1 (x) ωhx0 − xi ( x˙i0 − x˙i )(xl0 − xl ) dx0 dx. + ωhˇx − xi σi j hˇx − xi δ p j − dˇx Mlp + 2 ω0 (x) ω0 (ˇx) |ˇx − x|2

=

x∈B x0 ∈B

xˇ ∈B

(51)

Comparing Eqs. (17) and (51), the force state is obtained as "Z #  1   ηj ξpξ j  ωhηi 1 1 −1 ti (x)hηi = + (x) ωhηi ηl . σ hηi + ωhξi σi j hξi δ p j − 2 dξ Mlp ω0 (x) i j |η|2 2 H(x) ω0 (x) ω0 (x + ξ) |ξ|

(52)

The improvements in model stability using the proposed bond-associated approach over the initial state-based model will be shown in the numerical examples (see Section 7.1). To ensure the balance of angular momentum be maintained, the net moment exerted by t at x is evaluated using Eq. (52), the symmetry in σ, and ω0 (x) = ω0 : Z Mi = i jk η j tk dη, ZH Z Z  η j ηp ξn ξ p  ωhξi ωhηi −1 = dξ, i jk σkp hηi 2 dη + i jk ωhηi η j ηm dη Mmn σkp hξi δnp − ω0 |η| |ξ|2 H H ω0 |H {z } Z

Z

M jm

 η j ηp ξn ξ p  ωhξi ωhηi dξ, σkp hηi 2 dη + i jk δ jn σkp hξi δnp − ω0 |η| |ξ|2 H H ω0 Z Z   η j ηp η j ηp ωhηi ωhηi = σkp hηi 2 dη + σkp hηi δ jp − dη, i jk i jk ω ω |η| |η|2 0 0 H H Z ωhηi = i jk σ jk hηi dη, H ω0 = 0, =

i jk

thus the balance of angular momentum holds for any influence state that only depends on relative distance between material points (not limited to spherical functions). 210

6. Correspondence damage modeling at the bond level Tupek et al. (2013) presented a state-based damage correspondence model for peridynamics. In the state-based failure model, a material point loses all its associated bonds once its damage reaches a critical value. Such approach can lead to issues in numerical simulations due to loss of a finite amount of material. As an alternative method, we propose a bond-associated damage correspondence theory such that the failure criterion is enforced at the bond-level. Therefore, in a damaged area only bonds that meet the damage criterion are broken (or degraded), rather than removing material points. In theory, not only does such a bond-associated damage modeling approach improve stability, but also it results in localized damage (see Sections 7.4 and 7.5 for numerical demonstrations). Inspired by the state-based damage framework (Tupek et al., 2013, Section 2), the influence state is modified as   ˆ |η|, qhηi , ωhηi = ω where qhηi is the set of bond-associated internal variables, e.g. plastic strain, temperature, porosity, stress triaxiality, Lode angle, etc. and is governed through classical constitutive relations. A simple, practical assumption is to describe the material damage by a single scalar parameter D. The influence function is then formed as  

ω ˆ |η|, qhηi = ωη |η| ωD Dhηi , 13

215

220


where D ∈ [0, 1] is the bond-associated damage parameter, and ωη |η| is the conventional, spherical influence function which describes the relative degree of interactions between two material points associated with an undamaged bond. To incorporate damage irreversibility, ωD is required to be a non-increasing function of Dhηi. An additional requirement is that ωD (0) = 1 and ωD (1) = 0. These conditions insure that the degree of interactions of the associated material points are intact for the undamaged material, decreases as damage accumulates, and vanishes completely when the bond is fully damaged. A gradual decrease in the influence state for the transition zone can be prescribed. Alternatively, an immediate bond breakage approach from an undamaged to damaged bond can be taken. Definition 1. For post-processing/visualization purposes, damage at a material point x is defined as the degraded weighted volume of its associated bonds, given by R

ω |η| ωD Dhηi dη H(x) η R Φ(x) = 1 − . (53)
ω |η| dη H(x) η

A material point with damage Φ = 0 at a given time has intact interaction levels with all the material points that position in its attributed kernel space. A fully damage material point with Φ = 1, on the other hand, does not interact with any of its neighboring points. The specific form we adopt for the spherical influence function is   !2    |η|  1      if 0 < |η| < δ,
 2 1 − δ  |η| ωη |η| =      0 otherwise.

The damage-related part of the influence function is applied through a continuum failure model as discussed in the sequel. 6.1. Brittle failure

225

Brittle fracture occurs in small deformations, with plasticity only occuring in a small region near the crack tip satisfying the assumptions of linear elastic fracture mechanics and the notion of small scale yielding. It typically involves immediate crack growth and unstable crack propigation. To agree with the nature of fast rupture in brittle materials, it is appropriate to bypass the gradual transition in influence function degradation and break a bond at once when it meets a given failure criterion, e.g. a Tresca criterion, a widely-used classical models for describing brittle failure. Let’s propose the following bond damage parameter    ˆ ≤ σcr , 0 if σhηi Dhηi =   1 otherwise,

where σ ˆ is the maximum equivalent stress (can be Tresca, von Mises, or any other form of equivalent stress) the bond η has experienced over the course of loading, and σcr is the fracture strength. 6.2. Ductile fracture Failure in ductile materials often involves a gradual process in material degradation. To allow a transitional decay in the loadcarrying capacity of a bond, the following bond damage parameter is defined, which involves a linear decrease in the influence function   1 − Dhηi   if Dhηi > Dc ,
  1 − Dc ωD Dhηi =     1 otherwise,

where Dc ≤ 1 is a material constant (Dc = 1 corresponds to a sudden bond breakage). Alternative forms of decay, e.g. quadratic, can also be introduced. Plastic strain and stress triaxiality are often noted as key factors in deriving ductile failure. In addition, loading rate can play an important role depending on the nature of problem. The Johnson-Cook failure criterion (Johnson and Cook, 1985) is a well-accepted classical damage theory for describing ductile fracture. A Johnson-Cook correspondence model is employed to evolve the bond damage parameter as follows Z t

Dhηi =

˙ Dhηi dτ,

0

 p ε˙ hηi      ˙ εf Dhηi =    0

if Dhηi < 1, otherwise,

14

where ε p hηi is the bond-associated equivalent plastic strain, and ε f is defined as !# " " # σm hηi ε˙ p 1 + d4 ln ε f = d1 + d2 exp −d3 , σe hηi ε˙ 0

230

where σm hηi and σe hηi are the bond-associated hydrostatic and von Mises equivalent stresses, respectively, and d1 , d2 , d3 , d4 , and ε˙ 0 are material constants. Effects of other parameters, e.g. temperature, can be incorporated in a similar manner. 7. Numerical examples The following numerical examples are provided for multiple purposes: 1.) to observe the enhancements in stability achieved by the means of the bond-associated approach compared to the state-based model, 2.) to compare simulation results with analytical solutions or numerical results obtained from other mathematical theories, and 3.) to showcase the capabilities of the proposed model in problems involving damage. The semi-Lagrangian theory is implemented in Peridigm (Parks et al., 2012), an open-source peridynamic code, and simulations were run on the Texas Advanced Computing Center (TACC) Stampede supercomputer. CUBIT mesh generation software (Blacker et al., 1994) is used to create geometries and their descretized forms. The introduced bond-associated approach is used in all the simulations, and the numerical algorithm of Flanagan and Taylor (1987) is employed for integration of the constitutive equations in rate form using a co-rotational stress rate to ensure material objectivity. Elastic response is modeled using a simple linear response, and J2 plasticity governs the plastic behavior. To prevent a strong shock loading, essential boundary conditions are applied through a prescribed velocity condition that is stationary initially and increases smoothly, i.e.,    π t     if t < t0 , v 1 − cos 0  2 t0 (54) v(t) =    v0 if t ≥ t0 .

235

240

245

For visualization purposes, the nodal values, e.g. stress and strain tensors, are evaluated using the state-based notation, i.e. updating nodal values using Eq. (20). Note that the nodal values do not play a role in solving the peridynamic equation of motion, in which only the bond-associated quantities contribute.

Remark 10. In theory, there is a bond between each material point. It can be very costly (or practically impossible in large simulations) to consider every possible bonds in numerical simulations. To reduce the computational cost, a search bucket (circle in 2D and sphere in 3D) can be used with a radius greater than the horizon to initialize bonds. In the following examples, we choose the search radius to be about 1.5-2 times the horizon. The bond-associated properties are tracked in a total Lagrangian way throughout a simulation. Note that the influence state is still determined with respect to the current positions. Therefore, the bonds outside an updated kernel are not involved in computing the nodal velocity gradient tensor. In problems with extremely large distortions, e.g. shock wave propagation, a much larger search radius is needed. While this approach increases the required memory allocation and increases computational cost, it eliminates potential complications associated with updating the neighbors for a peridynamic node, e.g. how to transfer the bond-associated properties from the previous set of neighbors to the new ones. An alternative approach might be to use a smaller search radius and update the neighborhoods periodically. However, a careful study on how to redistribute bond-associated properties (e.g. strain, stress, and damage), with an emphasis on consistency with the laws of thermodynamics, is required first. This reflection will be answered in future work. 7.1. On the material stability

250

255

A simple experiment is designed to highlight the effect of the material instability inherent in the state-based model on the simulated behavior of material. Consider a 3D peridynamic bar under uniaxial tensile loading where velocity boundary conditions are prescribed to its ends. The peridynamic equation of motion is solved using an explicit time integration scheme, once using the state-based model and another time using the bond-associated approach. As shown in 3, unstable modes of deformation are evident in the state-based results, leading to large errors in the predicted bulk response of the material. The interested reader can refer to the stability analysis of Behzadinasab and Foster (2020) with more discussions on the effects of a similar material instability. 7.2. Plate with a hole problem The rectangular plate with a small hole, Kirsch’s problem (Kirsch, 1898), is considered here. Suppose an elastic 2D infinite plate with a central hole of radius a is subjected to a far-field uniaxial loading, as shown in Fig. 4. Analytical solution for this classical case exists as follows: " #  a 2  a 2  a 4 ! σ ˜ σrr (r, θ) = 1− + 1−4 +3 cos2θ . 2 r r r 15

(a) Macroscopic response

(b) State-based displacements

(c) Bond-level displacements

Figure 3: Simulation results of a uniaxial test (along the y-axis) of a 3D semi-Lagrangian cuboid. The predicted macroscopic behavior (a) indicates an unstable behavior using the state-based model. The instabilities appearing in the displacements (mid-place y-values) simulated using the state-based model (b) are corrected using the bond-associated approach, resulting in a smooth gradient (c).

" #  a 2  a 4 ! σ ˜ 1+ − 1+3 cos2θ . 2 r r  a 2  a 4 ! σ ˜ −3 σrθ (r, θ) = − 1 + 2 sin2θ. 2 r r

σθθ (r, θ) =

Figure 4: Schematic of infinite plate with a small hole, subjected to a far-field uniaxial loading.

260

265

As it is impossible to simulate the fully infinite scenario, the problem is approximated with a finite plate with much larger dimensions than the hole radius; in this case, a square plate with an edge 12 times bigger than the hole diameter. Since peridynamics does not directly deal with traction boundary conditions, the far-field conditions are applied as a body force density on a volumetric region near the ends (width of the region is set equal to the horizon). Explicit dynamic time integration is employed until a set amount of time, then the stress values are normalized by the applied far-field value. The peridynamic simulations were run at two different node densities, labeled as coarse and fine in the following figures, with 10 and 20 nodes along the hole diameter, respectively. The horizon is set to 0.3 times the diameter (3 times the coarse nodal spacing). Figs. 5 and 6 provide a comparison between the simulation results and the exact solution. In the figures, dimensions are normalized by the hole radius. As shown in Fig. 5, σrθ contours are in good agreement between the analytical solution and the numerical model (the refined case). The peridynamic results of σrr (r, θ = 0) also compare well with the exact values, as shown 16

in Fig. 6. Small differences are expected to be related to the approximation of an infinite plate with a finite one.

(a) Simulation (refined discretization)

(b) Exact

Figure 5: Good agreement between the simulated and analytical contours of σrθ .

Figure 6: Simulated σrr along the x-axis is in good agreement with the analytical solution.

7.3. Taylor impact test 270

275

280

The Taylor impact experiment (Taylor, 1948) has been broadly used as a benchmark problem for validating numerical techniques involving large plastic deformations. The test involves a high-velocity impact of a cylindrical metallic bar with a rigid wall. A schematic of the experiment is shown in Fig. 7. Previous peridynamic studies (Foster et al., 2010; Tupek et al., 2013) modeled the experiment within the framework of the Lagrangian state-based peridynamics without damage. In order to compare results with the peridynamic simulations of Tupek et al. (2013), in addition to other numerical results from the literature (Kamoulakos, 1990; Zhu and Cescotto, 1995; Camacho and Ortiz, 1997; Li et al., 2010), the specific configuration given in (Simo, 1992) is utilized. The considered cylinder is 32.4 mm long with a radius of 3.2 mm. It is made of copper, which is modeled as a linear isotropic hardening material with properties provided in Table 1. The impact velocity is 227 m/s. The short-range repulsive forces approach of (Silling and Askari, 2005) is employed to model the contact between the cylinder and the wall. The computations are carried out using an explicit dynamic solver. To see the effects of numerical discretization in the calculations, three different cases of 8540 (D1), 18144 (D2), and 35956 (D3) peridynamic nodes are considered. 17

Figure 7: Schematic of the Taylor impact experiment. Table 1: Material properties used to simulate the Taylor impact test.

Properties Undeformed density Young’s modulus Poisson’s ratio Initial yield stress Linear hardening modulus

285

Values 8930 kg/m3 117 GPa 0.35 400 MPa 100 MPa

Fig. 8 depicts the distribution of equivalent plastic strains at the final configurations of the simulated bar. Table 2 collects the final length of the bar (l f ), the final radius of the impacting face of the bar (rmax ), and the maximum equivalent plastic strain p (max ) predicted by the semi-Lagrangian peridynamics and compares with the results of other researchers. We observe that our results compare reasonably well with those obtained using the finite element analysis (Kamoulakos, 1990; Zhu and Cescotto, 1995; Camacho and Ortiz, 1997) (FEM calculations performed in axisymmetric mode), the optimal-transportation meshfree approximation (Li et al., 2010), and the Lagrangian state-based peridynamics (Tupek et al., 2013).

(a) D1

(b) D2

(c) D3

Figure 8: Equivalent plastic strain contours at the simulated final configurations of the Taylor impact test.

7.4. Mode-I crack growth This example demonstrates modeling of brittle crack growth using the proposed framework. As shown in Fig. 9, a 2D plate with an initial crack is exposed to prescribed velocities over its left-side bounding region. Note that all the bonds that 18

Table 2: Taylor impact test of copper specimen. Comparison of simulation results.

Kamoulakos (1990) Zhu and Cescotto (1995) Camacho and Ortiz (1997) Li et al. (2010) Tupek et al. (2013) Semi-Lagrangian PD, D1 Semi-Lagrangian PD, D2 Semi-Lagrangian PD, D3

290

l f (mm) 21.47–21.66 21.26–21.49 21.42–21.44 21.43 21.4–21.5 21.27 21.37 21.40

R f (mm) 7.02–7.12 6.89–7.18 7.21–7.24 6.8 7.1–7.5 6.29 6.53 6.99

p max 2.47–3.24 2.75–3.03 2.97–3.25 3.0 2.69–3.29 2.25 2.84 3.25

initially cross the crack region are explicitly broken at the beginning. The Tresca failure model is adopted as the criterion for bond breakage. A smooth displacement field and an expected damage pattern are observed using an explicit time integration as displayed in Fig. 10.

Figure 9: Initial configuration for the mode-I crack growth example. Upward and downward velocities are prescribed over the red and blue regions, respectively.

(a) Displacement

(b) Damage

Figure 10: Mode-I crack growth. Displacements are magnified by a factor of 20.

7.5. Ductile failure This last example serves as an illustrative example for ductile fracture modeling using the semi-Lagrangian framework. A 3D specimen with dimensions 150 mm x 50 mm x 5 mm includes a hole with radius 10 mm and is discretized with a grid spacing 19

of 1 mm. The horizon size is chosen to be 3 mm. The bar is modeled using a correspondence isotropic hardening material with Young’s modulus 195 GPa and Poisson’s ratio 0.27. The flow stress σy is defined using the following powerlaw hardening rule:  ε p n σy = σ0 1 + ε0

295

300

where ε p is the equivalent plastic strain, and the material constants are set to σ0 = 311 MPa, ε0 = 0.0002, and n = 0.178. The Johnson-Cook correspondence failure criterion introduced in Section 6.2 is used to model failure in this ductile material. The damage model constants are set to d1 = 0.25, d2 = 0.3, and d3 = −1.5. A uniaxial tensile test is simulated by applying velocity boundary conditions at the ends of the material. A dynamic time integration scheme is used in the simulation. Fig. 11 shows the equivalent plastic strain contours at the crack initiation and the damage values at the complete failure. The predicted macroscopic response of the material is shown in Fig. 12. The simulations resemble a typical set of results for such an experiment. We want to emphasize that adopting a bond-associated damage modeling, rather than a state-based approach, resulted in a damage field localized to the most critical region of the material under loading. Using a state-based damage model would have resulted in a thick damage pattern (characterized by the size of the horizon).

(a) Crack initiation

(b) Complete failure

Figure 11: Simulated equivalent plastic strain at crack initiation (a) and damage at complete failure (b) in the ductile failure problem.

Figure 12: Macroscopic load-displacement response of the material in the ductile failure problem.

8. Conclusions and discussion

305

310

A semi-Lagrangian framework for peridynamics, which defines the material point interactions based on their positions in the current deformed configuration, has been presented. In contrast to the proposed model here, the Lagrangian peridynamics uses the initial undeformed configuration to define the kernel function, making it unsuitable for modeling extreme events, e.g. shock waves and earth moving simulations, in which significantly large material deformation occurs. Evaluation of bond forces in the proposed model is achieved in a convenient way through computation of a nonlocal velocity gradient L¯ – an intermediate quantity – then incorporating well-established local constitutive equations. Stability analysis indicates that a pure state-based view involves material instabilities. A bond-level formulation has been introduced to enhance the stability. Material failure is 20

315

320

325

330

achieved through a bond-level correspondence damage modeling approach. Numerical examples have been provided to illustrate the capabilities of the new theory. Besides simulation of material failure, contact modeling is a primary concern for many practitioners of peridynamics, e.g. in impact and penetration problems. Modeling of contact in the peridynamic theory has been mostly addressed through the short-range forces approach of Silling and Askari (2005), in which contact forces are analogous to repulsive spring forces and two additional material constants are involved. The introduced semi-Lagrangian framework opens up new pathways to contact modeling in peridynamics. For example, a natural contact model has been developed for semi-Lagrangian RKPM (Chi et al., 2015), which evaluates contact forces very similarly to internal forces. The natural contact formulation uses the same constitutive relations for material and contact models and thus does not involve any additional contact model parameters. Similarly, a peridynamic natural contact model can be developed within the semi-Lagrangian formulation. A nonlocal frictional contact model was also recently presented for Lagrangian peridynamics (Kamensky et al., 2019), based upon an approximate velocity gradient tensor. The friction model can be readily adopted within the semi-Lagrangian framework. Further research is required to more deeply probe the robustness of the proposed semi-Lagrangian formulation in modeling extreme events, in which the Lagrangian theory becomes unreliable. A Lagrangian finite deformation peridynamic model was recently utilized, in the context of Sandia Fracture Challenge 2017 (Kramer et al., 2019), to predict ductile fracture in additively manufactured metal (Behzadinasab and Foster, 2019) where large material deformations involved. While the approach resulted in a qualitatively good blind predictions, it quantitatively failed to accurately simulate the load-carrying capacity of the structure and its fracture response. It was reported that the material instability issues underlying the peridynamic model, as well as the Lagrangian nature of the peridynamic framework could have been the primary source of error in predicting the failure behavior of material. A revisit of these predictions with the presented semi-Lagrangian theory will be the subject of future work. Acknowledgments

335

The authors would like to acknowledge helpful discussions with Xiao Xu during the preparation of this manuscript. The authors are grateful for the financial support provided by the AFOSR MURI Center for Materials Failure Prediction through Peridynamics: Project NO. ONRBAA12-020. M. Behzadinasab also acknowledges the fellowship funding by The University of Texas at Austin. References References Amani, J., Oterkus, E., Areias, P., Zi, G., Nguyen-Thoi, T., Rabczuk, T., 2016. A non-ordinary state-based peridynamics formulation for thermoplastic fracture. Int. J. Impact Eng. 87, 83–94.

340

Behzadinasab, M., Foster, J. T., 2019. The third Sandia Fracture Challenge: peridynamic blind prediction of ductile fracture characterization in additively manufactured metal. Int. J. Fract. 218, 97–109. Behzadinasab, M., Foster, J. T., 2020. On the stability of the generalized, finite deformation correspondence model of peridynamics. Int. J. Solids. Struct. 182, 64–76.

345

Behzadinasab, M., Vogler, T. J., Peterson, A. M., Rahman, R., Foster, J. T., 2018. Peridynamics modeling of a shock wave perturbation decay experiment in granular materials with intra-granular fracture. J. Dynamic Behavior Mater. 4 (4), 529–542. Bergel, G. L., Li, S., 2016. The total and updated lagrangian formulations of state-based peridynamics. Comput. Mech. 58 (2), 351–370. Blacker, T. D., Bohnhoff, W. J., Edwards, T. L., 1994. Cubit mesh generation environment. volume 1: Users manual. Tech. rep., Sandia National Labs., Albuquerque, NM (United States).

350

Breitenfeld, M., Geubelle, P., Weckner, O., Silling, S., 2014. Non-ordinary state-based peridynamic analysis of stationary crack problems. Comput. Methods Appl. Mech. Eng. 272, 233–250. Breitzman, T., Dayal, K., 2018. Bond-level deformation gradients and energy averaging in peridynamics. J. Mech. Phys. Solids 110, 192–204.

355

Camacho, G., Ortiz, M., 1997. Adaptive lagrangian modelling of ballistic penetration of metallic targets. Comput. Methods Appl. Mech. Eng. 142 (3-4), 269–301. Chen, H., 2018. Bond-associated deformation gradients for peridynamic correspondence model. Mech. Res. Commun. 90, 34–41.

21

Chen, J.-S., Zhang, X., Belytschko, T., 2004. An implicit gradient model by a reproducing kernel strain regularization in strain localization problems. Comput. Methods Appl. Mech. Eng. 193 (27-29), 2827–2844.

360

Chi, S.-W., Lee, C.-H., Chen, J.-S., Guan, P.-C., 2015. A level set enhanced natural kernel contact algorithm for impact and penetration modeling. Int. J. Numer. Methods Eng. 102 (3-4), 839–866. Chowdhury, S. R., Roy, P., Roy, D., Reddy, J., 2019. A modified peridynamics correspondence principle: Removal of zero-energy deformation and other implications. Comput. Methods Appl. Mech. Eng. 346, 530–549. Flanagan, D., Taylor, L., 1987. An accurate numerical algorithm for stress integration with finite rotations. Comput. Methods Appl. Mech. Eng. 62 (3), 305–320.

365

Foster, J. T., Silling, S. A., Chen, W., 2011. An energy based failure criterion for use with peridynamic states. Int. J. Multiscale. Com. 9 (6). Foster, J. T., Silling, S. A., Chen, W. W., 2010. Viscoplasticity using peridynamics. Int. J. Numer. Methods Eng. 81 (10), 1242– 1258.

370

Foster, J. T., Xu, X., 2018. A generalized, ordinary, finite deformation constitutive correspondence model for peridynamics. Int. J. Solids. Struct. 141, 245–253. Gerstle, W., Sau, N., Silling, S., 2007. Peridynamic modeling of concrete structures. Nucl. Eng. Des. 237 (12-13), 1250–1258. Gingold, R. A., Monaghan, J. J., 1977. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon. Not. R. Astron. Soc. 181 (3), 375–389.

375

Ha, Y. D., Bobaru, F., 2010. Studies of dynamic crack propagation and crack branching with peridynamics. Int. J. Fract. 162 (1-2), 229–244. Jafarzadeh, S., Chen, Z., Bobaru, F., 2018. Peridynamic modeling of intergranular corrosion damage. J. Electrochem. Soc. 165 (7), C362–C374. Johnson, G. R., Cook, W. H., 1985. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures. Eng. Fract. Mech. 21 (1), 31–48.

380

Kamensky, D., Behzadinasab, M., Foster, J. T., Bazilevs, Y., 2019. Peridynamic modeling of frictional contact. J. Peridy. Nonloc. Model. 1 (2), 107–121. Kamoulakos, A., 1990. A simple benchmark for impact. Bench Mark, 31–35. Kilic, B., Madenci, E., 2010. Coupling of peridynamic theory and the finite element method. J. Mech. Mater. Struct. 5 (5), 707–733.

385

Kirsch, C., 1898. Die theorie der elastizitat und die bedurfnisse der festigkeitslehre. Zeitschrift des Vereines Deutscher Ingenieure 42, 797–807. Kramer, S. L., Boyce, B. L., Jones, A., et al., 2019. The third Sandia Fracture Challenge: predictions of ductile fracture in additively manufactured meta. Int. J. Fract. 218, 5–61.

390

Li, B., Habbal, F., Ortiz, M., 2010. Optimal transportation meshfree approximation schemes for fluid and plastic flows. Int. J. Numer. Methods Eng. 83 (12), 1541–1579. Littlewood, D. J., 2010. 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, pp. 209– 217.

395

Littlewood, D. J., 2011. A nonlocal approach to modeling crack nucleation in AA 7075-T651. In: ASME 2011 International Mechanical Engineering Congress and Exposition. American Society of Mechanical Engineers, pp. 567–576. Liu, W. K., Jun, S., Zhang, Y. F., 1995. Reproducing kernel particle methods. Int. J. Numer. Meth. Fl. 20 (8-9), 1081–1106. Madenci, E., Barut, A., Futch, M., 2016. Peridynamic differential operator and its applications. Comput. Methods Appl. Mech. Eng. 304, 408–451.

400

Mehrmashhadi, J., Chen, Z., Zhao, J., Bobaru, F., 2019. A stochastically homogenized peridynamic model for intraply fracture in fiber-reinforced composites. Compos. Sci Technol. 182, 107770. 22

Mitchell, J. A., 2011. A nonlocal ordinary state-based plasticity model for peridynamics. Tech. rep., Sandia National Laboratories (SNL-NM), Albuquerque, NM (United States). Mitchell, J. A., Silling, S. A., Littlewood, D. J., 2015. A position-aware linear solid constitutive model for peridynamics. J. Mech. Mater. Struct. 10 (5), 539–557. 405

Moës, N., Dolbow, J., Belytschko, T., 1999. A finite element method for crack growth without remeshing. Int. J. Numer. Methods Eng. 46 (1), 131–150. Nikravesh, S., Gerstle, W., 2018. Improved state-based peridynamic lattice model including elasticity, plasticity and damage. CMES-Comp. Modl. Eng. 116 (3), 323–347.

410

Oterkus, S., Madenci, E., Oterkus, E., Hwang, Y., Bae, J., Han, S., 2014. Hygro-thermo-mechanical analysis and failure prediction in electronic packages by using peridynamics. In: 2014 IEEE 64th Electronic Components and Technology Conference (ECTC). IEEE, pp. 973–982. Parks, M. L., Littlewood, D. J., Mitchell, J. A., Silling, S. A., 2012. Peridigm users’ guide v1. 0.0. Sandia report 7800. Seleson, P., Ha, Y. D., Beneddine, S., 2015. Concurrent coupling of bond-based peridynamics and the navier equation of classical elasticity by blending. Int. J. Multiscale Com. 13 (2).

415

Shojaei, A., Mudric, T., Zaccariotto, M., Galvanetto, U., 2016. A coupled meshless finite point/peridynamic method for 2d dynamic fracture analysis. Int. J. Mech. Sci. 119, 419–431. Silling, S. A., 2000. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids 48 (1), 175–209.

420

Silling, S. A., 2017. Stability of peridynamic correspondence material models and their particle discretizations. Comput. Methods Appl. Mech. Eng. 322, 42–57. Silling, S. A., Askari, E., 2005. A meshfree method based on the peridynamic model of solid mechanics. Comput. Struct. 83 (17), 1526–1535. Silling, S. A., Bobaru, F., 2005. Peridynamic modeling of membranes and fibers. Int. J. Nonlin. Mech. 40 (2-3), 395–409.

425

Silling, S. A., Epton, M. A., Weckner, O., Xu, J., Askari, E., 2007. Peridynamic states and constitutive modeling. J. Elasticity 88 (2), 151–184. Silling, S. A., Lehoucq, R., 2010. Peridynamic theory of solid mechanics. In: Advances in applied mechanics. Vol. 44. Elsevier, pp. 73–168. Silling, S. A., Parks, M. L., Kamm, J. R., Weckner, O., Rassaian, M., 2017. Modeling shockwaves and impact phenomena with eulerian peridynamics. Int. J. Impact Eng. 107, 47–57.

430

Simo, J. C., 1992. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Comput. Methods Appl. Mech. Eng. 99 (1), 61–112. Swegle, J., Hicks, D., Attaway, S., 1995. Smoothed particle hydrodynamics stability analysis. J. Comput. Phys. 116 (1), 123–134. Taylor, G. I., 1948. The use of flat-ended projectiles for determining dynamic yield stress I. theoretical considerations. Proc. R. Soc. Lond. A 194 (1038), 289–299.

435

Tupek, M. R., Radovitzky, R., 2014. An extended constitutive correspondence formulation of peridynamics based on nonlinear bond-strain measures. J. Mech. Phys. Solids 65, 82–92. Tupek, M. R., Rimoli, J. J., Radovitzky, R., 2013. An approach for incorporating classical continuum damage models in statebased peridynamics. Comput. Methods Appl. Mech. Eng. 263, 20–26.

440

Warren, T. L., Silling, S. A., Askari, A., Weckner, O., Epton, M. A., Xu, J., 2009. A non-ordinary state-based peridynamic method to model solid material deformation and fracture. Int. J. Solids. Struct. 46 (5), 1186–1195. Yaghoobi, A., Chorzepa, M. G., 2015. Meshless modeling framework for fiber reinforced concrete structures. Comput Struct 161, 43–54. Zaccariotto, M., Mudric, T., Tomasi, D., Shojaei, A., Galvanetto, U., 2018. Coupling of FEM meshes with peridynamic grids. Comput. Methods Appl. Mech. Eng. 330, 471–497.

445

Zhu, Y., Cescotto, S., 1995. Unified and mixed formulation of the 4-node quadrilateral elements by assumed strain method: application to thermomechanical problems. Int. J. Numer. Methods Eng. 38 (4), 685–716. 23