Concurrent coupling of peridynamics and classical elasticity for elastodynamic problems

Accepted Manuscript Concurrent coupling of peridynamics and classical elasticity for elastodynamic problems Xiaonan Wang, Shank S. Kulkarni, Alireza Tabarraei

PII: DOI: Reference:

S0045-7825(18)30470-5 https://doi.org/10.1016/j.cma.2018.09.019 CMA 12079

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date : 12 April 2018 Revised date : 25 July 2018 Accepted date : 13 September 2018 Please cite this article as: X. Wang, et al., Concurrent coupling of peridynamics and classical elasticity for elastodynamic problems, Comput. Methods Appl. Mech. Engrg. (2018), https://doi.org/10.1016/j.cma.2018.09.019 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

Highlights (for review)

•

A coupled peridynamic-finite element method is developed to model dynamic fracture.

•

Spurious wave reflections is effectively suppressed.

•

Arlequin method is used to impose the mechanical compatibility between finite element and peridynamic sub regions.

*Manuscript Click here to download Manuscript: Peridynamics1.pdf

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Click here to view linked Referenc

Concurrent coupling of peridynamics and classical elasticity for elastodynamic problems Xiaonan Wanga , Shank S. Kulkarnia , Alireza Tabarraeia,∗ a Department of Mechanical Engineering and Engineering Science, University of North Carolina at Charlotte, Charlotte, NC 28223, USA

Abstract In this paper, a coupled peridynamic-finite element method for modeling mechanical behavior and damage growth in materials is developed. Peridynamics is a nonlocal continuum model which, in contrast to other continuum models, uses integration instead of spatial derivations in its governing equations. Utilizing integration instead of derivatives is advantageous in modeling fracture since the governing equations remain valid even after the initiation or growth of discontinuities. Although peridynamics can capture material fracture effectively, however due to its nonlocal formulation peridynamics is computationally expensive. To reduce the computational costs, we propose to couple peridynamics with finite elements and use peridynamics only in small zones where higher accuracy is needed. The main challenge in developing such a coupling method is to eliminate the artifacts introduced by the interface of the two subdomains. One of the main issues is spurious wave reflections which occurs because high frequency waves traveling from peridynamic region cannot enter the finite element zone and spuriously reflect back into the peridynamic zone. This will lead to an increase in the energy of the peridynamic zone and will drastically reduce the computational accuracy. The main feature of the proposed method is eliminating the spurious wave reflections such that the coupled method is as accurate as the pure peridynamics. Keywords: Peridynamics, Finite Elements, Arlequin method, Dynamic fracture

1. Introduction A problem of significant importance in solid mechanics is the accurate modeling of fracture. Several numerical techniques have been developed in the past for treating discontinuities and predicting crack growth. Element deletion method is one of the simplest finite element (FE) techniques for fracture 5

modeling. While the method is called element deletion method, the element is not deleted, but when the fracture criteria are satisfied in an element, the stress in the element is set to zero, i.e. zero ∗ Corresponding

author Email address: [email protected] (Alireza Tabarraei)

Preprint submitted to Elsevier

July 25, 2018

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

material resistance [1, 2]. Although its simplicity has made it one of the most widely used technique in industry, however this method performs poorly in predicting crack speed and growth path [2]. Nodal release techniques [3, 4, 5] and nodal splitting techniques [6, 7, 8, 9] are two other early finite 10

elements techniques developed for modeling crack propagation when the crack growth path is known a priori. In these techniques the domain is discretized such that the crack path coincides with the element boundaries. Crack growth in an arbitrary direction would require remeshing to realign element boundaries with the direction of crack growth. Remeshing is a cumbersome process and transfer of data from the original mesh to the new mesh results in the loss of accuracy.

15

Extended finite element method (X–FEM) [10] has been developed to alleviate the shortcomings associated with remeshing. In the X–FEM, enrichment functions are added to the finite element approximation using the framework of partition of unity [11]. This enables the domain to be modeled without explicitly meshing the crack surfaces, hence allowing arbitrary crack growth without remeshing. While extended finite element method has been used extensively in modeling both quasi-static [12,

20

13, 14, 15, 16] and dynamic fracture [17, 18, 19], tracking the crack surfaces within the framework of X–FEM is challenging which limits the applicability of X–FEM for 3D crack modeling. Another class of numerical methods used for fracture modeling are meshless methods (MMs). Meshless methods have been developed primarily to eliminate the issues associated with the presence of a mesh. Opposed to techniques such as finite element or finite volume, in meshless methods the

25

approximation is built from nodes only. Several meshless techniques including smooth particle hydrodynamics (SPH) [20, 21, 22], element–free Galerkin (EFG) [23] and reproducing kernel particle [24] have been developed. The main disadvantages of meshless methods are the higher computational costs and the treatment of essential boundary conditions. Since the shape functions of meshless methods are rational, accurate integration requires higher order integration schemes. Also because the shape

30

functions do not satisfy the Kronecker delta property, the imposition of essential boundary conditions is not as straightforward as in finite elements. More recently, fracture modeling techniques based on peridynamics (PD), a nonlocal continuum theory have been developed [25, 26]. In peridynamics, the forces acting at a particle are obtained by an integral operator that sums pairwise internal forces exerted on the particle by all the particles located

35

at a finite distance from the particle of interest. In contrast to classical continuum theory, peridynamics theory does not make any assumption on the differentiability of displacement field for obtaining the forces. The main advantage of peridynamics is that no assumption is made on the continuity of the displacement field hence discontinuity in the displacement field due to the presence of cracks does not necessitate special treatments. Peridynamics performance has been validated by applying it to several

40

sophisticated applications including polycrystals fracture [27], fracture of composite materials [28, 29, 30, 31, 32], structure stability and failure analysis [33], modeling of structure response under extreme 2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

loading [34], material fragmentation under impact [35] dynamic fracture analysis [36, 37, 38, 39, 40], simulation of the kinetic of phase transformation [41] and modeling heat transition in bodies with evolving cracks [42]. In addition to the numerical verifications, rigorous mathematical analysis have 45

been used to examine the properties of peridynamics. Silling and Lehoucq [43] showed that in the limit when nonlocal region around a point goes to zero, peridynamics converges to classical elasticity theory. Du and Zhou [44] developed a functional analytical framework for peridynamics and demonstrated the connections between peridynamics and classical elastic theory. Alali and Lipton [45] analytically investigated the multiscale dynamics of heterogeneous media using peridyanmic formulation.

50

Although the peridynamic theory is capable of modeling damage formation and growth without resorting to any external criteria, numerical simulations using peridynamics are computationally expensive. The high computational cost is attributed to the nonlocality of the theory. In nonlocal theories each particle interacts with a large group of particles, resulting in costly assembly operations of the nonlocal discrete systems. To reduce the computational consts, Lindsay et. al [46] proposed

55

decomposing the domain into overlapping subdomains and using different time steps in different subdomains. Smaller time step is used where higher accuracy is required and larger time steps are used elsewhere. Another issue associated with peridynamics is the prescription of displacement and traction boundary conditions. Since the variation formulation of the peridynamics does not include tractions, the forces acting on the surface are prescribed as body forces acting within a fictitious boundary layer

60

under the surface. Similar to tractions, displacements boundary conditions are imposed by constraining the displacement of material points within a fictitious boundary layer. A linear interpolation is used to approximate the value of the displacement of material points in the boundary layer based on the boundary conditions and the displacement of the points within the domain. Therefore the imposition of boundary conditions in peridynamics is not as accurate as techniques such as finite element method.

65

In this paper, to alleviate the disadvantages associated with peridynamics, we propose a technique for coupling peridynamics with finite elements (FE) for dynamic fracture modeling. In our proposed coupling technique, peridynamics description is used only in the areas of the domain where nucleation or growth of discontinuities is probable, and finite elements are used elsewhere. Therefore, fine scale behavior is captured by the peridynamics zone and finite elements are used in the zones where solution

70

is smooth. The proposed method takes advantages of the salient features of both techniques i.e. straight forward application of boundary conditions and lower computational cost of finite element modeling along with the superiority of peridynamics in modeling cracks growth and nucleation. The main challenge associated with developing coupled peridynamics–finite element methods is to minimize the fictitious interface effects. In static problems, special considerations need to be taken

75

to remove ghost forces, while in dynamic problems, the spurious wave reflections at the interface of peridynamics and finite element domains must also be eliminated. The main source of wave reflections 3

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

can be attributed to the change of constitutive behavior between an inherently nonlocal peridynamic domain and a local continuum domain, and the significant difference between the resolution of finite element mesh and the resolution of peridynamic zone. Since the minimum wave length which can 80

be supported by FE is much larger than the minimum wave length which can be resolved by the PD zone, the FE/PD interface acts as a rigid boundary for those components of the wave which can not be resolved by the finite element mesh. So instead of passing smoothly into the FE zone, the short wavelength (high frequency) components of the wave will reflect back into the PD domain. This will lead to spurious growth of the energy of the peridynamic domain, and will drastically reduce the

85

computational accuracy. The problem of spurious wave reflection at the interface of peridynamics and finite elements is studied analytically in Reference [47]. Coupling between peridynamics and finite element has received attention recently. Macek and Silling [48] proposed a technique for implementing the PD formulation in a conventional FEM code. The used truss elements to peridynamic bonds and used a “fuzzy” zone to impose the displacement

90

constraint between the FEM and PD zones. Kilic and Madenci [49] proposed a method for coupling peridynamics with finite element for elastostatic problems by employing an overlapping zone between peridynamic and finite element zone. Liu and Hong [50] used interface elements to couple finite element and peridynamic zones for static analysis. Lee et. al [51] used this approach to implement a parallel peridynamic–finite element code. Lubineau et. al [52] have used a morphing strategy to couple a

95

peridynamic zone to a classical continuum domain for static problems. Erkan et. al [53] studied a submodeling approach in which the boundary conditions of PD are imported from FEM . However they assume that the submodeling details do not affect the FEM simulation and the boundary of the PD is far enough from the local features. Galvento et. al [54, 55] have developed a coupling method between peridynamic and finite element zone by modifying the stiffness matrix of the finite element

100

zone. Li et. al [56] proposed a method to couple bond based peridynamics with FEM using node sharing technique. In this approach, both finite element and peridynamic zones should be discretized using the same size node spacing which imposes a restriction on the finite element mesh. In this paper, we propose a technique for coupling peridynamic and finite element subregions. In developing our coupling techniques we aim at satisfying the following two main criteria:

105

1. The peridynamic zone should be properly glued to the finite element subdomain such that a smooth enough wave can travel between subdomains without getting distorted at the interface. 2. Spurious wave reflections should be eliminated. For this purpose, the high frequency waves (fine scale oscillations) which cannot travel into the finite element zone should be appropriately damped.

110

In this paper, we propose a concurrent coupling technique in the framework of Arlequin approach 4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[57, 58], which uses Lagrange multipliers to glue two disparate subdomains. The Arlequin approach has been used previously for coupling atomistic zone to continuum zone [59, 60, 61, 62, 63] and to couple a nonlocal continuum domain to a local continuum zone [64]. In this paper, this approach is used to couple two continuum zones to each other. The main advantage of this proposed method is it 115

can be used to static as well as dynamic problems. Also in this technique, the discretization sizes of peridynamic and finite element zones are independent from each other, hence it is not required that the finite element mesh be refined to the PD node spacing at the interface of the two zones. To simplify the formulations, we present the coupling for linear elastic materials, however the proposed method is generic and is not restricted to linear elastic materials.

120

2. Formulation of elastodynamics in classical continuum Consider an elastic body occupying domain Ωc ⊂ Rd , where d is the number of space dimensions. Let Γ denotes the boundary of Ωc . The boundary Γ is decomposed of two nonoverlapping subregions denoted with Γt and Γu such that Γ = Γt ∪ Γu .

(1)

The strong form of the initial/boundary value problem for elastodynamics is as follow ∇ · σ + ρb = ρ¨ u in Ωc ,

(2)

where ∇ is the gradient operator, σ is the stress tensor, ρ is the mass density, b is the body force per ¨ signifies two time derivatives of the displacement field. unit mass, u is the displacement field and u We consider the following initial conditions for the displacement and velocity fields u(x, 0) = u0 (x) x ∈ Ωc ,

(3a)

˙ u(x, 0) = v0 (x) x ∈ Ωc ,

(3b)

where x denotes the Euclidean coordinates. Furthermore, we assume the following boundary conditions

¯ on Γu , u=u

(4a)

n · σ = ¯t on Γt ,

(4b)

¯ and ¯t are the the prescribed displacement and where n is normal vector to the boundary Γ and u traction vector field on portion Γu and Γt of the boundary respectively. The weak form is obtained by using the principle of virtual work by assuming the domain undergoes an arbitrary displacement field δu which gives rise to compatible virtual strains δ and can be written as

Z

Ωc

T

¨ dΩ + ρδu u

Z

δ : σdΩ =

Ωc

Z

Γt

5

T

δu tdΓ +

Z

Ωc

δuT bdΩ.

(5)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The weak form is discretized in space using finite elements. Assuming the finite element mesh consists of NF nodes, the displacement and strains are approximated as uh (x) =

NF X

NF X

NI (x)uI and δuh (x) =

I=1

h (x) =

NF X

NF X

BI (x)uI and δh (x) =

I=1

125

NI (x)δuI ,

(6a)

BI (x)δuI ,

(6b)

I=1

I=1

where uI is the displacement vector of node I, NI and BI are respectively the shape function matrix and strain–displacement matrix of node I. By using Eq. (6) in Eq. (5) and exploiting the arbitrariness of the virtual displacement field we obtain ¨ I = fIext mcI u




c

− fIint




c

,

(7)

where mcI is the lumped mass of node I and the vectors of external forces (fIext )c and internal forces
fIint c of node I are given by Z Z
fIext c = NI ρbdΩ + NI tdΓ, (8a) Ωc Γt Z
fIint c = BTI σdΩ. (8b) Ωc

2.1. Hamiltonian of classical elasticity

The total Hamiltonian of a classical continuum domain is due to the kinetic energy Wckin , internal strain energy Wcint and external potential energy Wcext and is given by Hc = Wckin + Wcint − Wcext .

(9)

The kinetic energy is given by Wckin =

NF X pcI · pcI , 2mcI

(10)

I=1

where pcI = mcI u˙ I is the momentum of node I and u˙ is the velocity vector. The internal energy is Z 1 Wcint = σij ij dΩ. (11) 2 Ωc The external potential energy is given by Wcext =

NF X I=1

fIext




c

uI .

(12)

The vectors of internal and external forces acting at node I can be obtained using the internal and external energies using


∂Wcint , ∂uI
∂Wcext fIext c = . ∂uI fIint

c

=

6

(13a) (13b)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 1: A peridynamic domain before and after deformation.

3. Basic formulation of peridynamics An alternative approach to classical continuum mechanics is peridynamic theory [25, 35]. Peridynamic is a nonlocal continuum theory in which a point x interacts with other points in its vicinity. As shown in the Fig. 1, the interaction zone of point x is its spherical neighborhood Hx = {x0 | 0 < kx − x0 k ≤ δ}, where δ is the horizon of point x. The vector from a point x to point x0 ∈ Hx

is called a bond defined by ξ = x0 − x and the bond length is denoted by ξ = |ξ|. The deformation

state of a bond ξ at time t is defined as Y[x, t] hξi = y(x + ξ, t) − y(x, t),

(14)

y(x, t) = x + d(x, t),

(15)

where

and d is the displacement vector field. Y[x, t] acts on bond ξ and produces the image of the bond under deformation. The extension scalar state e measures the change in the bond length due to deformation

e = y − x = |η| ,

y = |y0 − y| = |ξ + η| ,

The volume dilatation θ at point x is obtained from Z 3 3 θ(e) = (ωx) • e = ωxydVx0 − 3, q q Hx

x = |ξ| .

(16)

(17)

where q is the weighted volume defined by q = (ωx) • x = 7

Z

Hx

ωx2 dVx0 ,

(18)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where Hx denotes the horizon of point x and ω is the influence function which depends only on |ξ|, 130

ω = 1 for an unbroken bond and ω = 0 for a broken bond. The • operator denotes the inner product between two states and is defined in [26]. The deviatoric part ed of extension state e is defined by ed = e −

θx . 3

(19)

3.1. Equations of motion and numerical discretization The peridynamic equation of motion of point x is Z ¨ ρ(x)d(x, t) = (T[x, t] hx0 − xi − T[x0 , t] hx − x0 i)dVx0 + b(x, t) = fpint (x, t) + b(x, t),

(20)

Hx

where T is the force state. The force state T[(x, t)] takes the bonds connected to point x as input and produces a force density vector as the output. Since T acts on a vector and produce a vector, T is 135

similar to a tensor with the difference that T is not necessarily linear or continuous. fpint (x, t) is the equivalent of the classical internal force acting at point x. Numerical implementation of peridynamic continuum model requires the discretization of the domain. The most common discretization scheme employed in peridynamics is the meshfree method proposed in [35]. Opposed to the finite element discretization, in meshefree method the domain is

140

discretized by nodes instead of elements. In a meshfree discretization nodes are not connected to each other by elements or any other geometrical constraints. The discretized form of the peridynamic equations of motion of node I is

¨ I , t) = ρ(xI )d(x

L X

(T[xI , t] hxJ − xI i − T[xJ , t] hxI − xJ i)VcJ + b(xI , t) = fIint

J=1




p

(t) + b(xI , t), (21)

where xI denotes a peridynamic discrete node, xJ is a node whose volume overlaps with the horizon of xI and L is the total number of nodes whose volume overlaps with the horizon of xI . VcJ is the portion of volume of xJ which in the initial configuration is within the horizon of xI , as shown in Fig. 2a. VcJ can be given as VcJ = vIJ VJ , where VJ is the volume of xJ in the initial configuration and vIJ is the volume fraction given by    1   δ−|xJ −xI | 1 vIJ = 2 + rp     0

if

|xJ − xI | ≤ δ −

if

δ−

if

rp 2

rp 2 ,

< |xJ − xI | < δ +

|xJ − xI | ≥ δ +

rp 2 ,

(22)

rp 2

where rp represents the spacing between nodes. As shown in Fig. 2, volume fraction vIJ continuously

reduces from 1 to 0. The discretized equations of motion can be integrated in time using an explicit 145

time integration scheme such as the velocity Verlet method. 8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 2: (a) Discretization of domain in to PD points (only half volume of point xJ is inside horizon of point xI ), (b) vIJ as a function of distance.

3.2. Linear elasticity using peridynamic formulation For elastic materials T only depends on Y([x, t]) and can be obtained from the strain energy density function ψ P D (Y) given by T(Y) = ∇ψ P D (Y),

(23)

where ∇ denotes the Fr´echet derivative with respect to Y. The strain energy density function of a linear elastic material is ψ P D (θ, ed ) =

Kθ2 15µ + (ωed ) • ed , 2 2q

(24)

where K and µ are respectively the bulk modulus and shear modulus of material. Elastic materials can be represented in peridynamics using an ordinary model. In an ordinary model, T = tM, where M is the unit vector along the deformed bond direction and t is the magnitude of T. Using Eq. (24) in Eq. (23) the magnitude of the force state vector acting along the deformed bond direction is t=

3Kθ 15µ d ωx + ωe . q q

(25)

3.2.1. Plane stress peridynamic The volume dilatation in a plane stress state is [65] θ=

2(2v − 1) ωx • e , (v − 1) q

(26)

and the force state is [65] t=

2(2ν − 1) (ν − 1)



k0 θ −

β (ωed ) • x 3

9



ωx + βωed , q

(27)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

where β=

8µ , q

(28a) µ (v + 1)2 . 9 (2v − 1)2

0

k =k+

(28b)

3.3. Hamiltonian of peridynamics Assume the domain is discretized by N peridynamic nodes. Similar to classical continuum mechanics, the Hamiltonian of peridynamics is given by Hp = Wpkin + Wpint − Wpext ,

(29)

where the kinetic energy is defined by Wpkin

=

N X 1 I=1

2

d˙ TI ρd˙ I VI =

N X pp · pp I

I=1

2mpI

I

.

(30)

In this equation, d˙ I is the velocity of node I, VI is the volume of node I, mpI = ρVI and ppI = mpI d˙ I are respectively the mass and momentum of node I, and N is the total number of nodes in the peridynamic domain. The internal strain energy is N X

Wpint =

VI ψIP D ,

(31)

dI · bVI ,

(32)

I=1

and the potential energy of the external forces is Wpext =

N X I=1

150

where dI = d(xI , t) is the displacement of node I.
The force fIint p can be obtained from Wpint by [66] fIint




p

=

L X ∂Wpint = (T[xI , t] hxJ − xI i − T[xJ , t] hxI − xJ i)VcJ . ∂dI

(33)

J=1

Similarly, the external force acting at node I can be obtained from the potential energy of external forces using fIext




p

=

∂Wpext = bVI . ∂dI

(34)

4. Formulation of the Arlequin method In this section we describe a method based on the Arlequin framework [57, 58] for coupling finite elements with peridynamics. Following the Arlequin approach, the domain is subdivided into three subdomains: a pure peridynamics zone denoted by Ωp , a pure finite element zone (classical continuum 10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 3: Domain decomposition for Arlequin based coupling method.

155

zone) denoted by Ωc and an overlapping zone denoted by Ωo . The domain decomposition is shown in Fig. 3. To reduce the computational costs, the peridynamics is used only around cracks and finite elements is used elsewhere. In the Arlequin method the energy of the system is obtained as a linear combination of the energy from both models, i.e peridynamic elasticity and classical elasticity as ˙ = α(x)Hc (u, u) ˙ = H(u, u, ˙ d, d) ˙ + (1 − α(x))Hp (d, d)

α(x) Wckin + Wcint − Wcext + (1 − α(x)) Wpkin + Wpint − Wpext ,

(35a)

where u denotes the displacement field of the finite element zone and d represents the displacement field of the peridynamics zone. The weighting coefficients α should be chosen such that   1 ∀x ∈ Ω \Ω , c o α(x) =  0 ∀x ∈ Ω \Ω . p

(36)

o

The simplest form of coefficient α which satisfies the criteria of Eq. (36) is a piecewise linear function α(x) =

l1 l1 + l2

∀x ∈ Ωo ,

(37)

where l1 and l2 are respectively the distance from point x to the nearest points with α = 0 and α = 1.

The mechanical compatibility between the peridynamics and finite element zones requires that the classical continuum displacement u(x) be conforming with the peridynamics displacement field d(x) on the overlapping zone. The conformity between the displacements can be satisfied by defining constraints on the two displacement fields. Different form of constraints can be considered. For example, the following constraint requires that the displacement of peridynamic nodes should be equal to the displacement of finite elements at the location of peridynamics nodes gI = u(xI ) − dI =

X

J∈Mo

NJ (xI )uJ − dI = 0, 11

∀I ∈ No ,

(38)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 4: Lagrange multiplier interpolation.

160

where gI = {g1I , g2I , g2I } is the set of constraints imposed on displacement of node I, No is the set of all the peridynamic nodes in the overlapping zone and Mo is the set of all the finite element nodes whose support intersects the overlapping zone. This constraint is applied using Lagrange multipliers. For this purpose, the total Hamiltonian of the system is modified to HL (u, d, λ) = α(x)Hc (u) + (1 − α(x))Hp (d) +

X

I∈No

λI · g I ,

(39)

where λTI = {λ1I , λ2I , λ2I } is the Lagrange multipliers at peridynamic node I. Lagrange multipliers are approximated using a λ–mesh which is constructed over the overlapping zone λI (x) =

X

NJλ (xI )λJ ,

(40)

J∈Q

where NJλ is the shape function of node J of the Lagrange multiplier filed, Q is the set of Lagrange multiplier nodes and λJ is the value of Lagrange multiplier at λ–node J. Using Eq. (38) and Eq. (40) in Eq. (39) we obtain HL (u, d, λ) = α(x)Hc (u) + (1 − α(x))Hp (d) + Hλ (u, d, λ), where Hλ (u, d, λ) =

X X

λ NK (xI )λK

I∈No K∈Q

·

X

J∈S

NJ (xI )uJ − dI

!

(41)

.

(42)

A Lagrange multiplier mesh along with the finite element mesh and peridynamic grid is shown 165

in Fig. 4. Any mesh which satisfies the Ladyzhenskaya-Babuska-Brezzi (LBB) condition can be used to construct Lagrange multipliers shape functions. We consider two types of λ–mesh. A λ–mesh whose nodes coincides with the finite element nodes, i.e NJλ (xI ) = NJ (xI ), and a λ–mesh whose nodes coincide with the peridynamic nodes, i.e NJλ (xI ) = δIJ where δIJ is the Kronecker delta. We refer to the former method as a weak constraint coupling method and the second method a strict constraint

170

coupling method. 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

4.1. Equations of motion We use the Hamilton’s relations to find the equations of motion of finite element and peridynamic nodes. 4.1.1. Equations of motion of finite element nodes The Hamilton’s equations for a finite element node I in the classical continuum zone Ωc can be written as ∂HL , ∂uI ∂HL αI u˙ I = , ∂pcI αI p˙ cI = −

(43a) (43b)

where αI = α(xI ). By using Eq. (41) in Eq. (43) we obtain
∂ (α(x)Wcext ) ∂ α(x)Wcint ∂Hλ ∂(α(x)Hc ) ∂Hλ − = − − , =− ∂uI ∂uI ∂uI ∂uI ∂uI ∂Hc ∂Wckin αI u˙ I = αI = α . I ∂pcI ∂pcI

αI p˙ cI

By using Eq. (42) we obtain

X ∂Hλ = GIK λK , ∂uI

(44a) (44b)

(45)

K∈Q

where

X

GIK =

λ NK (xJ )NI (xJ ).

(46)

J∈No

By defining the scaled vector of external forces at finite element node I as Z Z
∂ (α(x)Wcext ) Fext = = α(x)N ρbdΩ + α(x)NI tdΓ, I I c ∂uI Ωc Γt
Z
∂ α(x)Wcint int FI c = = α(x)BTI σdΩ, ∂uI Ωc

(47a) (47b)

and using Eqs. (47) and (45) in Eq. (43) we obtain αI p˙ cI = Fext I αI u˙ I = αI




c

− Fint I

pcI . mcI




c

−

X

GIK λK ,

(48a)

K∈Q

(48b)

By combining Eq. (48) the equations of motion of finite element nodes are obtained as ¨ I = Fext MIc u I




c

− Fint I




c

−

X

GIK λK ,

(49)

K∈Q

where MIc is the scaled mass of finite element node I defined as MIc = αI mcI . 13

(50)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

175

4.1.2. Equations of motion of peridynamic nodes The equations of motion of peridynamic nodes can be obtained similarly. The Hamilton’s equation for a peridynamic node I can be written as ∂HL , ∂dI ∂HL (1 − αI )d˙ I = . ∂ppI

(1 − αI )p˙ pI = −

(51a) (51b)

Using Eq. (41) in Eq. (51) yields ∂Wpext ∂Wpint ∂Hp ∂Hλ ∂Hλ − = (1 − αI ) − (1 − αI ) − , ∂dI ∂dI ∂dI ∂dI ∂dI ∂Wpkin ∂Hp . (1 − αI )d˙ I = (1 − αI ) p = (1 − αI ) ∂pI ∂ppI

(1 − αI )p˙ pI = −(1 − αI )

(52a) (52b)

By using Eq. (42) we obtain ∂Hλ = −λI . ∂dI

(53)

Using Eqs. (33), (34) and (53) in Eq. (51) gives (1 − αI )p˙ pI = (1 − αI ) fIext (1 − αI )d˙ I = (1 −

pp αI ) Ip . mI




p

− (1 − αI ) fIint




p

− λI ,

(54a) (54b)

By defining the scaled peridynamic forces as
= (1 − αI ) fIext p ,

Fint = (1 − αI ) fIint p , I p

Fext I




p

(55a) (55b)

and combining Eq. (54), the equations of motion of peridynamic node I can be written as ¨ I = Fext MIp d I




p

− Fint I

where MIp is the scaled mass of peridynamic node I




p

MIp = (1 − αI )mpI .

− λI ,

(56)

(57)

Eqs. (49) and (56) provide the complete equations of motion for finite elements and peridynamic nodes. 4.2. Explicit Time Integration The time integration is conducted by using the velocity Verlet algorithm. Since the value of Lagrange multipliers at the beginning of each time step is not known we use a prediction–correction method similar to the technique developed in [60] to update the nodal displacements and velocities. 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In this algorithm the displacements at step n + 1 are obtained using the displacement, velocity and acceleration of step n by 1 n 2 ¨ ∆t un+1 = unI + u˙ nI ∆t + u I 2 I 1 ¨n 2 dn+1 = dnI + d˙ nI ∆t + d ∆t I 2 I

∀I ∈ M,

(58a)

∀I ∈ N .

(58b)

where M and N are set of all the finite element nodes and all the peridynamic nodes respectively. The half step velocities are calculated using n+ 21

u˙ I

1

n+ d˙ I 2

∆t n ¨ u 2 I ∆t ¨ n = d˙ nI + d 2 I = u˙ nI +

∀I ∈ M,

(59a)

∀I ∈ N .

(59b)

The trial accelerations at step n+1 are calculated by ignoring the force due to the Lagrange multipliers as i
1 h ext
n+1 int n+1 F − F I I c c MIc h i
n+1
1 int n+1 − F = p Fext I I p p MI

¨ ∗n+1 u = I

∀I ∈ M,

(60a)

¨ ∗n+1 d I

∀I ∈ N .

(60b)

The trial velocities at step n + 1 are obtained from n+ 12

u˙ ∗n+1 = u˙ I I

1

n+ d˙ ∗n+1 = d˙ I 2 I


∆t ∗n+1 ∆t ∗n ¨I ¨I + u ¨ ∗n+1 u = u˙ nI + u ∀I ∈ M, I 2 2   ∆t ¨ ∗n+1 ∆t ¨ ∗n ¨ ∗n+1 + dI = d˙ nI + dI + dI ∀I ∈ N . 2 2 +

(61a) (61b)

Since no Lagrangian multiplier force acts on the nodes not located in the overlapping zones, the 180

displacement and velocities of such nodes are obtained directly using Equations 60 and 61. The correct velocities of the nodes located within the overlapping zone are obtained by considering the Lagrange multipliers force in the calculation of accelerations as " # ∆t 1 X 1 X n n+1 n+1 n+1 n n ¨I − c ¨I − c ∀I ∈ Mo u˙ I = u˙ I + u GIK λL + u GIK λL 2 MI MI L∈Q L∈Q (62a) ∆t X n+ 12 ∗n+1 = u˙ I − c GIK λL ∀I ∈ Mo , MI L∈Q " # X X ∆t 1 1 n n+1 n+1 n+1 n n λ λ ¨ − p ¨ d˙ I = d˙ I + d Nk (xI )λL + d − Nk (xI )λL ∀I ∈ No I I 2 MI MIp L∈Q

L∈Q

(62b)

∆t X λ n+ 1 = d˙ ∗n+1 − p Nk (xI )λL 2 I MI L∈Q

∀I ∈ No ,

15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 5: The displacement of peridynamic nodes is discomposed into fine and coarse components. n+ 12

where λL g˙ In+1 =

=

X

J∈S

=

X

J∈S

1 2

 n  n+1 . By using Eq. (62) in the time derivative form of Eq. (38) we obtain λ L + λL

NJ (xI )u˙J n+1 − d˙ n+1 I NJ (xI )u˙J ∗n+1 −

∀I ∈ No

X

NJ (xI )

X

NJ (xI )

(63) ∆t X λ ∆t X n+ 12 n+ 12 ∗n+1 ˙ − d + , G λ N (x )λ JK L I L k I MJc MIp L∈Q

J∈S

L∈Q

which can be simplified as g˙ In+1 = g˙ I∗n+1 −

∆t X λ ∆t X n+ 12 n+ 1 λ + G NL (xI )λL 2 , JK L MJc MIp

J∈S

L∈Q

(64)

L∈Q

where g˙ I∗n+1 is defined as ˙ I )∗n+1 − d˙ ∗n+1 g˙ I∗n+1 = u(x = I

X

J∈S

NJ (xI )u˙ ∗n+1 − d˙ ∗n+1 J I

∀I ∈ No ,

(65)

Since g˙ In+1 = 0, Eq. (64) can be written in compact form as X

where AIK =

AIK λL = g˙ I∗n+1 ,

X ∆t ∆t NJ (xI )GJK − p NLλ (xI )I, c MJ MI

(66)

(67)

J∈S

where I is the identity matrix. To save in the computational costs, matrix A is diagonalized by " # X X ∆t ∆t λ AII = NJ (xI )GJK − p NL (xI ) . (68) MJc MI L∈Q

J∈S

The Lagrange multiplier of each λ-node is obtained from ˙ I∗n+1 . λI = A−1 II g

(69)

5. Removing the spurious reflections The method presented in Section 4 provides mechanical coupling between the classical elasticity and peridynamics zones and significantly reduces the spurious wave reflections. In this section, we improve the efficiency of the method in removing the spurious wave reflections by adding a damping 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

term to the equations of motion of the peridynamic nodes located within the overlapping zone. For this purpose, we decompose the oscillations of the peridynamic nodes into fine and coarse oscillations as shown in Fig. 5. The decomposition can be written as dI = dcoarse + dfine I I

∀I ∈ No .

(70)

The course scale oscillations can be resolved by the finite element mesh, but the fine scale oscillations cannot be resolved by the finite element mesh and will be spuriously reflected back by finite element mesh if they are not properly damped. For the purpose of consistency of deformation between peridynamic and finite element zone, we assume that the coarse scale displacement of each peridynamic node is equal to the displacement of the classical elasticity subdomain at the location of the peridynamic node. Using the finite element shape functions the coarse scale displacements can be approximated as dcoarse = I

X

NJ (xI )uJ

J∈S

∀I ∈ No .

(71)

By combining Eq. (71) and Eq. (70) the fine scale displacements are dfine = dI − I

X

NJ (xI )uJ

∀I ∈ No ,

(72)

NJ (xI )¨ uJ

∀I ∈ No .

(73)

J∈S

and subsequently the fine scale oscillations are ¨ fine = d ¨I − d I

X

J∈S

By multiplying both sides of Eq. (73) with MIp and using Eqs. (49) and (56) we obtain ¨ fine = Ffine MIp d I I

∀I ∈ No ,

(74)

where Ffine is I Ffine = Fext I I




p

X p

− MI

− Fint I NJ (xI )




p

+

X

λ NK (xI )λK

K∈Q

(Fext J )c

− Fint J

J∈S




− MJc c

P

K∈Q

GJK λK

(75) .

Fine scale oscillations cannot be resolved by the finite element mesh. Therefore, they can not transmit into the finite element zone and will reflect back if they are not damped appropriately. We damp the fine scale oscillations by modifying the equations of motion given by Eq. (74) to include a viscous damping term as ¨ fine = Ffine − M p CI d˙ fine MIp d I I I I

∀I ∈ No ,

(76)

where d˙ fine is the velocity of fine scale oscillations obtained from I d˙ fine = d˙ I − I

X

NJ (xI )u˙ J

J∈S

17

∀I ∈ No ,

(77)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

and CI is the damping coefficient. In this paper, we approximate the damping function using the following parabolic equation as CI = Cvw

αI2 , LB

(78)

where C is a user defined constant, vw is the wave speed at current time t, and LB is the overlapping length. Hence the value of viscous damping is zero at the interface of pure peridynamic and the 185

overlapping zone and gradually increases to Cvw at the interface of overlapping zone and pure finite element domain. A small value for C cause insufficient damping of fine oscillations while a large value of C leads to numerical instability [67]. In this paper, a value of C = 30 is adopted [63, 14] in solving the numerical examples. The velocity Verlet algorithm is used to integrate the fine scale equation of motion (Eq. (74)) in time  fine 1  fine
fine fine ¨n dn+1 = (dnI ) + ∆t d˙ nI + ∆t2 d , I I 2       fine fine fine 1 n+ 1 ¨n = d˙ nI + ∆t d , d˙ I 2 I 2    fine
1 n+ 21 fine n+1 fine ˙ ¨ n+1 F − C d d = , I I I I MIp fine 1   fine  fine 1 ˙ n+ 2 ¨ n+1 d˙ n+1 = d ∆t d . + I I I 2

(79a) (79b) (79c) (79d)

After calculating the fine scale displacements and velocities, the total displacement of peridynamic nodes located in the overlapping zone is obtained by adding course and fine scale oscillations dn+1 = dIn+1 I

190


fine

+

X

J∈S

NJ (xI )uJ

∀I ∈ No .

(80)

6. Numerical Results In this section we present numerical results to verify the capability of the proposed method in removing the spurious reflections and providing a coupling between the peridynamics and classical elasticity zone. 6.1. One-dimensional longitudinal wave propagation

195

In the first example, the propagation of a longitudinal wave from the peridynamic zone to the finite element zone is modeled. The two–dimensional domain of a rectangular plate with a length of L1 = 400 cm and width of w = 8 cm is considered as shown in Fig. 6. Central portion of length L2 = 200 cm is modeled using peridynamics.

18

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 6: Configuration of the 2D rectangular domain.

Figure 7: Snapshots of the central line of plate. The height represent the displacement in the x direction. (a) Initial waves and (b) Wave position obtained using pure peridynamics modeling.

The plate has a Young’s modulus of 69 GPa, Poisson’s ratio of 0.3 and a mass density of 2700 kg/m3 . We assume that the plate is in a plane stress state. The center of the domain is at the origin and the initial displacement of the wave in the x direction is prescribed by    h i 2π a1 + cos 2π|x| + 1 + 1 + a2 cos (−0.5w1 ≤ x ≤ 0.5w1 ) , w1 w2 (x+0.5w1 ) dx |t=0 =  0 (x < −0.5w or x > 0.5w ) , 1

(81)

1

where dx represents the x–component of the displacement field, a1 = 0.0003 m, a2 = 0.3 m, w1 = 48cm

200

and w2 = 3cm . The initial wave is shown in Fig. 7a and is obtained by superposing a low frequency wave with high frequency waves. As a reference solution we first model the wave propagation using peridynamics description in the entire domain. The peridynamic equations of motions are integrated using a time step of ∆t = 10−8 s. The initial wave configuration and the wave configuration at t = 2.4 × 10−4 s are shown in Fig. 7a and

205

Fig. 7b respectively. As shown, the initial wave splits into two waves and each wave travels toward an edge of the domain with an amplitude equal to the half of the amplitude of the initial wave. To show the issue of spurious wave reflection, we solve the wave propagation problem by coupling finite elements and peridynamics without considering the remedies proposed in this paper for suppressing the spurious wave reflections. This coupling scheme is shown in Fig. 8, where the overlapping

210

zone is decomposed into two sub domains FE2PD and PD2FE. The coefficient α = 1 on the entire FE2PD zone and α = 0 over the whole PD2FE zone. In the FE2PD zone, the PD nodes are considered 19

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 8: A bad coupling between peridynamic and finite element which causes spurious wave reflection.

Figure 9: Snapshots of waves at t = 2.4 × 10−4 s. (a) ∆t = 10−7 and (b) ∆t = 10−8

as slaves which follows the displacement of the FE nodes and in the FE2PD zone the FE nodes are treated as slaves which follow the displacement of PD nodes. Coupling approaches similar to this has been used previously to couple finite elements and peridynamics. To study the impact of time step 215

size on the problem, the problems is solved using two time steps of ∆t = 10−7 s and ∆t = 10−8 s. The snapshots of the wave propagation at t = 2.4 × 10−4 s are shown in Fig. 9a and Fig. 9b. These figures show that when this ”bad coupling” scheme is used the low frequency component of the wave can travel into the finite element zone but the high frequency waves cannot travel into the finite element zone and reflects back into the PD zone. By decreasing the time step, the issue of wave

220

reflection becomes worse. Finally, we verify the performance of the proposed coupling method in eliminating the spurious reflections by modeling the longitudinal wave propagation. The impact of the Lagrange multipliers discretization on the results are studied by using two λ–meshes. In the first λ–mesh which is used, every λ–node coincide with a PD node while the second λ–mesh which is employed every λ–node coincides

225

with the finite element mesh. We will refer to the first technique as a strict compatibility method and to the second technique as a weak compatibility method. Snapshots of the wave propagation at t = 2.4 × 10−4 s using ∆t = 10−7 s and ∆t = 10−8 s are shown in Fig. 10a–d. These figures show that the weak compatibility method is more efficient in removing spurious reflections than the strict compatibility method. Comparing the snapshots of Fig. 10a and Fig. 10b show that by reducing

230

the integration time step the performance of strict compatibility method in removing the spurious reflections reduces. On the other hand, snapshots of Fig. 10c and Fig. 10d show that the performance

20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 10: Snapshots of waves at t = 2.4 × 10−4 s obtained using strict (a − b) and weak (c − d) compatibility formulation.

of the weak compatibility method in eliminating spurious reflections is not affected by the integration time step. To further study the wave propagation problem, the total Hamiltonian, kinetic energy and internal energy of the peridynamic zone are plotted in Fig. 11. If no wave reflection occurs the energy 235

of the PD zone should go to zero. Nonzero energy of the PD zone indicates that spurious reflections have occurred and high frequency waves are trapped in the PD zone. The energy of the PD zone versus time for the bad coupling method is shown in Fig. 11a and Fig. 11b. These figures show that when bad coupling is used the energy of the PD zone does not approach zero. A nonzero energy of the PD zone is due to the spurious reflections shown in Fig. 7a and Fig. 7b. Comparing Fig. 11a with

240

Fig. 11b shows that spurious reflections are more pronounced when a smaller integration time step is used. The energy of the PD zone when the finite element and PD zone are coupled using strict compatibility method is shown in Fig. 11c and Fig. 11d. Plots of Fig. 11c show that when a larger time step is used, the energy of the PD domain reduces considerably indicating that spurious reflections

245

are eliminated significantly. As shown in Fig. 11d the energy of the PD zone is higher when a smaller time step is used indicating that the strict compatibility method is more efficient in removing spurious reflections when the time step is larger. The total Hamiltonian, kinetic energy and internal energy of the PD zone when coupling is conducted using weak compatibility method are shown in Fig. 11e and Fig. 11f. These results show that

250

the energy of the PD zone goes to zero which indicates that weak method can successfully damp the high frequency waves and prevent spurious wave reflections. Moreover, the time step does not impact the performance of the weak compatibility method and even when the time step is small the high frequency waves are damped. Since the results of this example show that the weak compatibility method

21

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 11: Strain energy, kinetic energy and total Hamiltonian of the pure peridynamic zone.

has a superior performance in eliminating spurious reflections, we use this technique to solve the next 255

numerical examples. 6.2. Crack branching under mode I loading In this section we use the proposed coupling method to simulate the crack branching in a rectangular plate subjected to a tensile loading as shown in Fig. 12. The plate has a width of 0.1 cm and height of 0.12 cm and includes an edge crack with a tip at the center of the plate. The center part of the plate

260

is modeled using peridynamic formulation while finite element method is used in the top and bottom parts of the plate. The finite element zones are discretized by a mesh of uniform size of 0.004 cm. To study the effect of the PD size on the results we conduct the simulations by discretizing the PD zone using a coarse grid with uniform node spacing of 0.0005 cm and a fine grid with uniform node spacing of 0.00025 cm. For each case the overlapping zone is the size of 5 layers of finite elements. The

22

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 12: A rectangular plate with a horizontal edge crack. L1 = 0.024 cm, L2 = 0.02 cm, L3 = 0.315 cm, W = 0.1 cm.

265

plate is subjected to a uniaxial tensile stress of 60 MPa. The material parameters are the same as the example in Section 6.1. An integration time step of 1 ps is used for the purpose of time integration of the equations of motion for both finite element and peridynamic zones. The crack propagation path obtained using the proposed coupling method at t = 0.57 µs is shown in Fig. 13a and Fig. 13b. Both figures show similar crack growth path. The crack initially propagates

270

along a straight path and branches after that. These results are in agreement with the results obtained previously using extended finite element method [68] and pure peridynamic method [69]. 6.3. Modeling of impact In the third example the impact of a bar on a plate is modeled using the proposed coupling method by employing a weak compatibility. This example demonstrates the capability of this coupling method

275

is passing wave from finite element zone to peridynamics zone for a general three–dimensional problem.

23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

Figure 13: The crack branching path. Mesh size of PD is (a) 8 times finer than the FEM mesh and (b) 16 times finer than the FEM mesh.

As shown in Fig. 14, the projectile is a bar of length 10 cm with a square cross section of size 1 cm×1 cm. The Young’s modulus and Poisson’s ratio of the bar are respectively 200 GPa and 0.3 and its mass density is 7823 kg/m3 . The initial velocity of the bar is 100 cm/s and the target plate dimensions are 1.4 × 1.4 × 0.2 cm3 . As shown in Fig. 14 the middle part of the bar with a length of 4 cm is modeled 280

using peridynamics and the rest of the domain is discretized using a finite element mesh. The finite element mesh is composed of 8-node brick elements of size 1 mm and the PD zone is discretized by converting each brick element into 8 PD nodes. For the purpose of solving this problem we have linked our peridynamic code to the commercial finite element code Velodyne [70]. The velocity contour plots of the bar after hitting the plate are shown in Fig. 15. After the impact, a

285

compressive wave travels from the right end of the bar to the left end and the velocity of the bar goes to zero. These results indicate that the proposed coupling is capable of modeling the transfer of wave from finite element to peridynamic and from peridynamic to finite element for a general three-dimensional problem.

24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 14: Impact of a bar and plate using coupled peridynamics and finite element method. Peridynamic zone overlaps with 8 finite elements.

-4

t = 1.005 x10 s

-4

t = 1.105 x10 s

-4

t = 1.030 x10 s

-4

t = 1.130 x10 s

-4

-4

t = 1.055 x10 s

t = 1.155 x10 s

-4

t = 1.180 x10 s

-4

t = 1.080 x10 s

Figure 15: Contour velocities of the bar after impacting the plate. The velocity reduces after the bar hits the plate. This reduction of velocity is then spread from the right to the left of the bar. Proposed coupling method can transfer the wave from finite element to peridynamics and vice versa.

25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The results are further verified by comparing the velocity of point A (shown in Fig. 14) obtained 290

from the proposed coupling method and finite element method. A plot of velocity of point A versus time is shown in Fig. 16 and the results obtained from the coupling method are in good agreement with those obtained from the finite element method.

Figure 16: The time history of the velocity of a point within the right overlapping zone. The result from proposed coupling method is very close to the pure FEM simulation.

6.4. Crack propagation under mixed mode loading In this section, we use the proposed coupling method to model the crack propagation in a plate under mixed-mode fracture loading. As shown in Fig. 17, the domain is 3D specimen with dimensions of 6.25 × 6.25 × 1.6 cm3 . As shown in the figure, peridynamic description is used in the vicinity of the crack tip while the rest of the domain is modeled using finite elements. We assume the Young’s modulus is 29 GPa and Poisson’s ratio is 0.29. The critical energy release rate is G0 =36 KJ/m2 . The mixed-mode loadings is applied by prescribing the crack tip asymptotic displacement field given by r 1+υ r θ θ θ θ app ux = [KIapp cos (κ − 1 + 2sin2 )] + KII sin (κ + 1 + 2cos2 ), (82a) E 2π 2 2 2 2 1+υ uy = E

r

θ θ θ θ r app [KIapp sin (κ + 1 − 2sin2 )] − KII cos (κ − 1 − 2sin2 ), 2π 2 2 2 2

(82b)

to the boundary of the domain. Since the domain boundary is composed of finite element nodes, 295

the imposition of these boundary conditions is exact. In Eq. (82), ux and uy are the displacement app components in the x and y directions, KIapp and KII are applied stress intensity factors for mode

I and II loadings, r and θ are the polar coordinates with respect to the crack tip. The effective p app app 2 stress intensity factor of the mixed mode loading is defined as Keff = (KIapp )2 + (KII ) and the 26

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 17: The domain and crack geometry used for modeling crack propagation under mixed mode loading. The vicinity of the crack tip is modeled using peridynamic formulation. The FEM mesh is shown in green and peridynamics nodes are shown in blue. K app

loading phase angle describing the ratio of the mode I and II loadings is defined as, φ = tan−1 KIapp . II

300

The external displacement boundary conditions are applied by increasing the effective stress intensity √ factor in increments of ∆Keff = 103 MPa m per second, which is equivalent to a low loading rate of approximately 0.5 cm/s. A time step of ∆t = 5.0 × 10−8 s is used for the purpose in time integration of equations of motion in the entire domain. The crack propagation path for different loading phase angles are shown in Fig. 18. As expected

305

under pure mode I loading, i.e when φ = 0 the crack grows on a straight path and kinks when the phase angle is not zero. In Fig. 19 the kink angle obtained from the proposed coupling method is compared with the angle obtained from maximum circumferential stress criterion which assumes that the cracks propagate along the direction with minimum strain energy. The graphs of Fig. 19 show that the kink angles obtained using the proposed method is in good agreement with the theoretical

310

predictions. The stress intensity factor at which the crack propagates corresponds to the critical stress intensity factor of the domain and can be used to find the critical energy release rate using Gcr =

KI2 K2 + II . E E

(83)

The critical energy release rate obtained from our results is shown in Fig. 20. The results are in good

27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 18: Crack propagation path under quasi-static mixed-mode loadings. The crack propagation path is shown in red.

Figure 19: The crack kink angles obtained from the proposed coupling method and maximum circumferential stress criterion.

28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Figure 20: The critical energy release rate obtained from the simulations using the coupling method compared with the critical energy release rate of the material.

agreement with the energy release rate of the material.

7. Conclusion In this paper, an efficient technique for coupling peridynamics and finite element zones was pro315

posed. The proposed techniques can effectively remove the spurious reflections at the interface of finite element and peridynamic zones. This coupling method is based on the Arlequin approach in which the two domains overlap partially. Lagrange multipliers are used to enforce the compatibility of deformation between peridynamic and finite element zones. We examined the proposed method by verifying the performance of the method using two different λ–mesh. One of the λ–meshes is match-

320

ing the peridynamic nodes while the other one is matching the finite element mesh. Our numerical simulations showed that the latter mesh can damp the high frequency oscillations more efficiently and prevents the spurious reflections more effectively. Finally, the performance of the method was verified by modeling wave propagation and crack growth under different conditions.

8. Acknowledgments 325

The authors are grateful to Charlotte Research Institute for partial financial support of this project. The authors would like to thank Dr. Cameron Bell and Mr. Andrey Andreyev from Corvid Technologies for their help and support in using Velodyne. Simulations were performed on the University Research Computing’s (URC) High Performance Computing (HPC) clusters at UNC Charlotte. The authors would also like to thank anonymous reviewers for their valuable comments. 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

330

References [1] B. Bourdin, G. A. Francfort, J.-J. Marigo, Numerical experiments in revisited brittle fracture, Journal of the Mechanics and Physics of Solids 48 (4) (2000) 797–826. [2] J.-H. Song, H. Wang, T. Belytschko, A comparative study on finite element methods for dynamic fracture, Computational Mechanics 42 (2) (2008) 239–250.

335

[3] S. Aoki, K. Kishimoto, M. Sakata, Finite element computation of dynamic stress intensity factor for a rapidly propagating crack using ˆj-integral, Computational mechanics 2 (1) (1987) 54–62. [4] A. Kobayashi, A. F. Emery, S. Mall, Dynamic-finite-element and dynamic-photoelastic analyses of two fracturing homalite-100 plates, Experimental Mechanics 16 (9) (1976) 321–328. [5] G. Rousselier, Numerical treatment of crack growth problems, in: Advances in Elasto-Plastic

340

Fracture Mechanics, Ispra, Italy, Apr. 1979,, 1980, pp. 165–189. [6] J. Bakuckas, A. Lau, T. Tan, J. Awerbuch, Computational methodology to predict damage growth in unidirectional composites-I. theoretical formulation and numerical implementation, Engineering fracture mechanics 52 (5) (1995) 937–951. [7] P. O. Bouchard, F. Bay, Y. Chastel, I. Tovena, Crack propagation modelling using an advanced

345

remeshing technique, Computer Methods in Applied Mechanics and Engineering 189 (3) (2000) 723–742. [8] A. R. Khoei, R. Yasbolaghi, S. O. R. Biabanaki, A polygonal finite element method for modeling crack propagation with minimum remeshing, International Journal of Fracture 194 (2) (2015) 123–148.

350

[9] A. Shahani, M. A. Fasakhodi, Finite element analysis of dynamic crack propagation using remeshing technique, Materials & Design 30 (4) (2009) 1032–1041. [10] N. Moes, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International journal for numerical methods in engineering 46 (1999) 131–150. [11] J. M. Melenk, I. Babuˇska, The partition of unity finite element method: basic theory and appli-

355

cations, Computer methods in applied mechanics and engineering 139 (1-4) (1996) 289–314. [12] N. Sukumar, J.-H. Pr´evost, Modeling quasi-static crack growth with the extended finite element method part i: Computer implementation, International journal of solids and structures 40 (26) (2003) 7513–7537.

30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[13] P. R. Budarapu, R. Gracie, S. P. Bordas, T. Rabczuk, An adaptive multiscale method for quasi360

static crack growth, Computational Mechanics 53 (6) (2014) 1129–1148. [14] A. Tabarraei, N. Sukumar, Extended finite element method on polygonal and quadtree meshes, Computer Methods in Applied Mechanics and Engineering 197 (5) (2008) 425–438. [15] M. Xu, A. Tabarraei, J. T. Paci, J. Oswald, T. Belytschko, A coupled quantum/continuum mechanics study of graphene fracture, International journal of fracture 173 (2) (2012) 163–173.

365

[16] C. Comi, S. Mariani, Extended finite element simulation of quasi-brittle fracture in functionally graded materials, Computer Methods in Applied Mechanics and Engineering 196 (41) (2007) 4013–4026. [17] G. Zi, H. Chen, J. Xu, T. Belytschko, The extended finite element method for dynamic fractures, Shock and Vibration 12 (1) (2005) 9–23.

370

[18] I. Nistor, O. Pantal´e, S. Caperaa, Numerical implementation of the extended finite element method for dynamic crack analysis, Advances in Engineering Software 39 (7) (2008) 573–587. [19] J.-H. Song, T. Menouillard, A. Tabarraei, Explicit dynamic finite element method for failure with smooth fracture energy dissipations, Mathematical Problems in Engineering 2013. [20] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The astronomical journal

375

82 (1977) 1013–1024. [21] L. D. Libersky, A. G. Petschek, Smooth particle hydrodynamics with strength of materials, Springer Berlin Heidelberg, Berlin, Heidelberg, 1991, pp. 248–257. [22] L. D. Libersky, P. W. Randles, T. C. Carney, D. L. Dickinson, Recent improvements in SPH modeling of hypervelocity impact, International Journal of Impact Engineering 20 (6) (1997)

380

525–532. [23] T. Belytschko, Y. Y. Lu, L. Gu, Element-free Galerkin methods, International journal for numerical methods in engineering 37 (2) (1994) 229–256. [24] W. K. Liu, S. Jun, Y. F. Zhang, Reproducing kernel particle methods, International journal for numerical methods in fluids 20 (8-9) (1995) 1081–1106.

385

[25] S. Silling, Reformulation of elasticity theory for discontinuities and long-range forces, Journal of the Mechanics and Physics of Solids 48 (2000) 175–209. [26] S. A. Silling, M. Epton, O. Weckner, J. Xu, E. Askari, Peridynamic states and constitutive modeling, Journal of Elasticity 88 (2) (2007) 151–184. 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[27] E. Askari, F. Bobaru, R. Lehoucq, M. Parks, S. Silling, O. Weckner, Peridynamics for multiscale 390

materials modeling, Journal of Physics: Conference Series 125 (2008) 012078. [28] B. Kilic, A. Agwai, E. Madenci, Peridynamic theory for progressive damage prediction in centercracked composite laminates, Composite Structures 90 (2009) 141–151. [29] W. Hu, Y. D. Ha, F. Bobaru, Modeling dynamic fracture and damage in a fiber-reinforced composite lamina with peridynamics, International Journal for Multiscale Computational Engineering

395

9 (6) (2011) 707–726. [30] W. Hu, Y. Ha, F. Bobaru, Peridynamic model for dynamic fracture in unidirectional fiberreinforced composites, Computer Methods in Applied Mechanics and Engineering 217-220 (2012) 247–261. [31] J. Xu, A. Askari, O. Weckner, S. Silling, Peridynamic analysis of impact damage in composite

400

laminates, Journal of Aerospace Engineering 21 (2008) 187–194. [32] A. Askari, J. Xu, S. Silling, peridynamic analysis of damage and failure in composites, in: 44TH AIAA Aerospace Sciences Meeting and Exhibit, AIAA2006-8, pp. 1–12. [33] B. Kilic, E. Madenci, Structural stability and failure analysis using peridynamic theory, International Journal of Non-Linear Mechanics 44 (2009) 845–854.

405

[34] P. Demmie, S. Silling, An approach to modeling extreme loading of structures using peridynamics, Journal of Mechanics of Materials and Structures 2 (2007) 1921–1945. [35] S. Silling, E. Askari, A meshfree method based on the peridynamic model of solid mechanics, Computers and Structures 83 (2005) 1526–1535. [36] A. Agwai, I. Guven, E. Madenci, Predicting crack initiation and propagation using xfem, czm

410

and peridynamics: A comparative study, in: Electronic Components and Technology Conference (ECTC), 2010 Proceedings 60th, pp. 1178–1185. [37] S. Silling, O. Weckner, E. Askari, F. Bobaru, Crack nucleation in a peridynamic solid, International Journal of Fracture 162 (2010) 219–227. [38] B. Kilic, E. Madenci, Prediction of crack paths in a quenched glass plate by using peridynamic

415

theory, International Journal of Fracture 156 (2009) 165–177. [39] Y. Ha, F. Bobaru, Characteristics of dynamic brittle fracture captured with peridynamics, Engineering Fracture Mechanics 78 (2011) 1156–1158.

32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[40] Y. Ha, F. Bobaru, Studies of dynamic crack propagation and crack branching with peridynamics, International Journal of Fracture 162 (2010) 229–244. 420

[41] K. Dayal, K. Bhattacharya, Kinetics of phase transformations in the peridynamic formulation of continuum mechanics, Journal of the Mechanics and Physics of Solids 54 (2006) 1811–1842. [42] F. Bobaru, M. Duangpanya, A peridynamic formulation for transient heat conduction in bodies with evolving discontinuities, Journal of Computational Physics 231 (2012) 2764–2785. [43] S. Silling, R. Lehoucq, Convergence of peridynamics to classical elasticity theory, Journal of

425

Elasticity 93 (2008) 13–37. [44] K. Zhou, Q. Du, Mathematical and numerical analysis of linear peridynamic models with nonlocal boundary conditions, SIAM Journal on Numerical Analysis 48 (2010) 1759–1780. [45] B. Alali, R. Lipton, Multiscale dynamics of heterogeneous media in the peridynamic formulation, Journal of Elasticity 106 (2012) 71–103.

430

[46] P. Lindsay, M. Parks, A. Prakash, Enabling fast, stable and accurate peridynamic computations using multi-time-step integration, Computer Methods in Applied Mechanics and Engineering 306 (2016) 382–405. [47] S. Kulkarni, A. Tabarraei, An analytical study of wave propagation in a peridynamic bar with nonuniform discretization, Engineering Fracture Mechanics 190 (2018) 347–366.

435

[48] R. W. Macek, S. A. Silling, Peridynamics via finite element analysis, Finite Elements in Analysis and Design 43 (15) (2007) 1169–1178. [49] B. Kilic, E. Madenci, Coupling of peridynamic theory and the finite element method, Journal of Mechanics of Materials and Structures 5 (2010) 707–733. [50] W. Liu, J.-W. Hong, A coupling approach of discretized peridynamics with finite element method,

440

Computer Methods in Applied Mechanics and Engineering 245-246 (2012) 163–175. [51] J. Lee, S. E. Oh, J.-W. Hong, Parallel programming of a peridynamics code coupled with finite element method, International Journal of Fracture 203 (1-2) (2017) 99–114. [52] G. Lubineau, Y. Azdoud, F. Han, C. Rey, A. Askari, A morphing strategy to couple non-local to local continuum mechanics, Journal of the Mechanics and Physics of Solids 60 (2012) 1088–1102.

445

[53] E. Oterkus, E. Madenci, O. Weckner, S. Silling, P. Bogert, A. Tessler, Combined finite element and peridynamic analyses for predicting failure in a stiffened composite curved panel with a central slot, Composite Structures 94 (3) (2012) 839–850. 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[54] U. Galvanetto, T. Mudric, A. Shojaei, M. Zaccariotto, An effective way to couple FEM meshes and peridynamics grids for the solution of static equilibrium problems, Mechanics Research Com450

munications 76 (2016) 41–47. [55] M. Zaccariotto, T. Mudric, D. Tomasi, A. Shojaei, U. Galvanetto, Coupling of fem meshes with peridynamic grids, Computer Methods in Applied Mechanics and Engineering 330 (2018) 471–497. [56] H. Li, H. Zhang, Y. Zheng, H. Ye, M. Lu, An implicit coupling finite element and peridynamic method for dynamic problems of solid mechanics with crack propagation, International Journal

455

of Applied Mechanics 1850037. [57] H. B. Dhia, Multiscale mechanical problems:

the Arlequin method, Comptes Rendus de

l’Academie des Sciences Series IIB Mechanics Physics Astronomy 12 (326) (1998) 899–904. [58] H. B. Dhia, G. Rateau, The Arlequin method as a flexible engineering desing tool, International Journal For Numerical Methods in Engineering 62 (2005) 1442–1462. 460

[59] A. Tabarraei, X. Wang, A. Sadeghirad, J. H. Song, An enhanced bridging domain method for linking atomistic and continuum domains, Finite Elements in Analysis and Design 92 (2014) 36–49. [60] S. Xiao, T. Belytschko, A bridging domain method for coupling continua with molecular dynamics, Computer methods in applied mechanics and engineering 193 (17) (2004) 1645–1669.

465

[61] S. Prudhomme, H. B. Dhia, P. T. Bauman, N. Elkhodja, J. T. Oden, Computational analysis of modeling error for the coupling of particle and continuum models by the arlequin method, Computer Methods in Applied Mechanics and Engineering 197 (41) (2008) 3399–3409. [62] S. Zhang, R. Khare, Q. Lu, T. Belytschko, A bridging domain and strain computation method for coupled atomistic-continuum modelling of solids, International Journal for Numerical Methods in

470

Engineering 70 (8) (2007) 913–933. [63] A. Sadeghirad, A. Tabarraei, A damping boundary condition for coupled atomistic–continuum simulations, Computational Mechanics 52 (3) (2013) 535–551. [64] F. Han, G. Lubineau, Coupling of nonlocal and local continuum models by the Arlequin approach, International Journal For Numerical Methods in Engineering 89 (2012) 671–685.

475

[65] Q. V. Le, W. K. Chan, J. Schwartz, A two-dimensionalordinary, state-based peridynamic model for linearly elastic solids, International Journal for Numerical Methods in Engineering 98 (8) (2014) 547–561. doi:10.1002/nme.4642. 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[66] E. Madenci, E. Oterkus, Peridynamic theory and its applications, Vol. 17, Springer, 2014. [67] A. C. To, S. Li, Perfectly matched multiscale simulations, Physical Review B 72 (3) (2005) 035414. 480

[68] T. Belytschko, H. Chen, J. Xu, G. Zi, Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment, International Journal for Numerical Method in Engineering 58 (1-2) (2003) 1873–1905. [69] F. Bobaru, G. Zhang, Why do cracks branch? a peridynamic investigation of dynamic brittle fracture, International Journal of Fracture 196 (1-2) (2015) 59–98.

485

[70] Velodyne. URL http://corvidvelodyne.com/

35