A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads

Engineering Fracture Mechanics xxx (2017) xxx–xxx

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads Xiaoping Zhou a,b,c,⇑, Yunteng Wang a,b,c, Yundong Shou a,b,c, Miaomiao Kou a a

State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing 400045, PR China School of Civil Engineering, Chongqing University, Chongqing 400045, PR China c Key Laboratory of New Technology for Construction of Cities in Mountain Area, Chongqing University, Chongqing 400045, PR China b

a r t i c l e

i n f o

Article history: Received 20 January 2017 Received in revised form 18 July 2017 Accepted 24 July 2017 Available online xxxx Keywords: Linear elastic model Conjugated bond-based peridynamics Crack propagation Dynamic loads The standard bond-based peridynamics

a b s t r a c t A novel conjugated bond linear elastic model is proposed and implemented into the bondbased peridynamic (BB-PD) framework. In this model, micro-elastic PD bond energy is not only related to the normal stretch of bonds, but also related to the rotation bond angles of a pair of conjugated bonds. Therefore, micro bond energy mechanism in this study is different from that in the classical continuum mechanism or the standard BB-PD. Only one micro-elastic constant in the standard BB-PD results in the limitation of the effective Poisson’s ratio for isotropic material. However, the novel conjugated bond linear elastic model incorporating two micro-elastic constants are proposed, which can overcome the limitation of Poisson’s ratio in the standard BB-PD. By comparing the strain energy in the proposed model with that in the classical elastic model, the corresponding micromacro parameter relationships can be established. In addition, energy-based bond rupture criteria are implemented in the proposed numerical model to simulate fracture problems under dynamic loads. In order to verify the ability and accuracy of the proposed numerical model to simulate fracture problems under dynamic loads, some numerical examples are investigated. The present numerical results are in good agreement with the previous experimental and numerical results. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction In computation mechanics, numerical simulation of fracture problems for both static cracks and dynamic cracks is still an active and persistent challenge. There are some different methods to model cracks in finite element and meshfree methods. For the finite element method, one of the simple and robust methods is interelement separation model where cracks are modeled along element interfaces in the mesh [1–4]. Another simple method to model cracks was developed by Remmers et al. [5] who introduced crack segments in finite elements. The embedded discontinuity model was also proposed by Armero and Garikipati [6], Belytschko et al. [7], Samaniego et al. [8] in the finite element mesh to simulate crack problems. Based on the ‘local’ partition of unity, the extended finite element method (XFEM), which is a very accurate method for crack problems, was developed by Belytschko and Black [9] and Moes et al. [10]. For the previous meshfree method, Rabczuk and Belytschko [11,12] developed a ‘cracking particle’ model (CPM) in meshfree method, where discontinuities are introduced at

⇑ Corresponding author at: State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing 400045, PR China. E-mail address: [email protected] (X. Zhou). http://dx.doi.org/10.1016/j.engfracmech.2017.07.031 0013-7944/Ó 2017 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

2

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Nomenclature cðnÞ cn D E E f fn ft F G Gc HX I K mjik s sn uðXÞ _ uðXÞ € ðXÞ u w wl wb wn wc WC W PD W PD l W PD b x x0 X X0 d

eij hjik0 hjik

Hjik

kt

l

v

n

g q v

w

micro-modulus constant in bond-based peridynamics normal stiffness of bonds ratio of amount of broken bonds to the total amount of interaction in one horizon Young’s modulus elastic Lagrange strain tensor pairwise force function in conjugated bond-based peridynamics normal bond force density in conjugated bond-based peridynamics tangent bond force density approximate deformation gradient shear modulus fracture surface energy horizon of a given material point X identity tensor bulk elastic modulus moment vector density due to a pair of conjugated bonds bond stretch normal bond stretch displacement vectors of the material point X velocity vectors of the material point X acceleration vectors of the material point X scalar-valued micro-potential function normal stretch strain energy micro-potential rotation Strain energy micro-potential due to a pair of conjugated bonds micro-potential energy function critical strain energy potential in the bond strain energy density in classical elastic mechanics macro-elastic strain energy density macro-elastic stretch strain energy density macro-elastic rotation strain energy density location of a material point in deformed configuration neighbor particle for a given point x in deformed configuration location of a material point in the reference configuration neighbor particle for a given point X in the reference configuration radius of one horizon strain tensor in the bond nij conjugated bond angle in the reference configuration conjugated bond angle in the deformed configuration rotation angle of conjugated bonds in the deformed configuration rotation resistance stiffness scalar factor representing the broken of bond Poisson’s ratio bond vector in the reference configuration relative displacement vector between two interacting material points material density unit orientation vector of the bond unit orientation vector of the conjugated bond

the particle positions. The visibility criterion or some modification of it were also proposed in the meshfree methods for cracking [13–16]. Ventura et al. [17] enriched the MLS base functions around the crack tip to model kinked and curve cracks. Although these numerical methods have been successfully applied to simulate the fracture problems, theories of these methods are in the framework of continuum mechanics, whose governing equations of motion are based on partial derivatives with respect to the spatial coordinates. The assumption of continuity in continuum mechanics is inherently insufficient for modeling cracks as partial derivatives are undefined along the crack faces where the displacement field is discontinuous. Therefore, some computation methods involving displacement gradients or higher-order spatial derivatives in a domain containing a crack must remove the discontinuous displacement field by redefining the discretized body so that the crack lies on the boundary or by using other techniques for evaluating the spatial derivatives on crack surfaces [9,18,19]. However, the peridynamic formulation eliminates the spatial derivatives altogether by solely depending on an integral formulation of the force acting at a continuum point, resulting in equilibrium equations that are valid everywhere in the body [20].

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

3

Peridynamic (PD) theory is a nonlocal formulation of continuum mechanics, which was first proposed by Silling [21] to solve discontinuous problems, such as cracks, and long-range forces. The motion equation of PD theory is in an integraldifferential form instead of the spatial derivative terms, which naturally supports the presence of discontinuities in the deformation field. Another advantage of the PD theory is that the resulting equations of motion are naturally discretized using particle-based method [22]. PD formulation [21,23–25] can be generally divided into two categories, i.e. bondbased peridynamics (BB-PD) which is comparatively simple but restricted to simulate materials with the fixed Poisson’s ratio, and the state-based peridynamics (SB-PD) which overcomes the limitation of Poisson’s ratio by considering the effects of deformations of other material points in one horizon on two interacting material points [24,25]. Moreover, advantages of SB-PD make it possible to incorporate generally constitutive models [26–28]. Since BB-PD is comparatively simple, BB-PD was developed and applied to investigate some solid fracture problems. Silling and Askari [23] proposed the BB-PD formulation to simulate the deformation of a bar. Bobaru et al. [29] discussed the numerical convergence in the BB-PD theory by simulating the deformation of a bar. The dynamic fracturing problems were simulated using the BB-PD model by Silling [30] and Weckner and Abeyaratne [31]. Ha and Bobaru [32] and Agwai et al. [33] applied the BB-PD model to investigate the propagation and branching of cracks under dynamic loads. The characteristics of dynamic brittle fracture captured with BB-PD model were also studied by Ha and Bobaru [34]. Huang et al. [35,36] extended the BB-PD theory to quantitatively investigate propagation of cracks subjected to both quasi-static and dynamic loads. However, BB-PD theory is based on the assumption of pairwise interactions of the same magnitude, which leads to being restricted to Poisson’s ratio of 1/4 for isotropic 3D problems. Recently, Liu and Hong [37] proposed a force compensation scheme adding additional forces on each pairwise force between particles to model materials with Poisson’s ratio other than 1/4. Prakash and Seidel [38] decomposed the stretch of a bond into two parts along the orthogonal directions to overcome the restriction of the fixed Poisson’s ratio. Moreover, a BB-PD code was developed by Liu and Hong [39] for 3D simulation for brittle and ductile solids. Ha et al. [40] and Lee et al. [41] simulated the crack propagation and coalescence in rock-like materials under the uniaxial compression using the coupling FEM and BB-PD method with the force compensation scheme. Additionally, a novel contact algorithm was implemented for the interaction of PD regime and finite elements for the effective impact fracture analysis by Lee et al. [42]. Although some researchers developed and successfully simulated fracturing problems using BB-PD, several aspects of the framework of BB-PD need to be further studied. In the present study, a novel conjugated bond linear elastic model is proposed in the framework of BB-PD to simulate dynamic fracturing problems in the light of Stillinger-Weber potential functions [43] and the corrected lattice bond cell [44]. In contrast to the standard BB-PD, the micro strain energy potential in the conjugated BB-PD is divided into two parts: the normal stretch strain energy potential and the rotation strain energy potential of a pair of the conjugated bonds. Two micro bond constants, i.e. normal stiffness of bonds and rotation resistance stiffness of a pair of conjugated bonds, are derived using the equality relationship between micro strain energy density in the conjugated BB-PD and macro strain energy density in the classical elastic mechanics. The conjugated BB-PD can be considered as a typical case of the ‘‘bondpair” material model in Refs. [24,25], where only a pair of bonds with isometric and collinear characteristics is considered in computation. The conjugated BB-PD model has the advantage over the standard BB-PD model, which overcomes the limitation of a fixed Poisson’s ratio. Moreover, an energy-based damage model is incorporated into the conjugated BB-PD to model the initiation and propagation of cracks. In order to verify the ability and accuracy of the conjugated BB-PD, an isotropic plate with a circular cutout under quasi-static conditions is firstly simulated as a benchmark example. Then, three numerical examples subjected to dynamic loads are modeled using the proposed numerical method. Next, two types of numerical convergence are verified. Meanwhile, effects of mechanical properties, PD particle arrangements and PD particle orientations on the dynamic crack initiation and propagation are also studied to illustrate the characteristics of the proposed numerical method. This paper is organized as follows: The brief description of BB-PD theory is presented in Section 2. The conjugated bond linear elastic model is derived in Section 3. A benchmark example is performed to verify the ability and accuracy of the proposed numerical method in Section 4. Three dynamic fracturing problems are analyzed in Section 5. The conclusions are drawn in Section 6.

2. Bond-based peridynamic framework PD theory [21,23–25] is a non-local continuum model, but it is different from the conventional non-local models by replacing the spatial differential equations with the integral operations. PD theory is strongly non-local theory as opposed to a weakly non-local theory or a strain gradient type non-local theory. The internal force density is represented as an integral of forces instead of a divergence of stress in classical theory or an integral of stress in a weak non-local theory. Furthermore, the primary kinematic field variable is displacement as opposed to a strain or strain gradient which removes all continuity requirements, which results in the advantages of addressing discontinuous problems. In the bond-based peridynamic (BB-PD) model, every material point X in the continuum is assumed to interact with other material points X0 within a surrounding finite distance d, which is called the horizon HX , as shown in Fig. 1. Consider two interacting material points X and X0 in the reference configuration separated by the relative distance n ¼ X0  X. In the deformed configuration, two material points x and x0 are obtained by x ¼ X þ u and x0 ¼ X0 þ u0 , respectively. Thus, the relative displacement between two interacting material points is given by g ¼ u0  u and the relative position in the deformed Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

4

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 1. Schematic diagram in a BB-PD model.

configuration is n þ g. Within a horizon of X, material point X is interacted with other material points X0 through PD bonds. According to Newton’s law, PD equation of motion at the point X at time t is given by

qu€ ðX; tÞ ¼

Z HX

f ðX; X0 ; u; u0 ; tÞdV X0 þ bðX; tÞ kX0  Xk 6 d

ð1Þ

in which, f is the pairwise force function which is defined to describe the interaction between material points X and X0 within the horizon of X, b denotes the body force density at material point X. The pairwise force function contains all the constitutive information of materials, but it eliminates any spatial derivatives or assumption on the continuity of displacement. The pairwise force function f can be reduced to

f ¼ f ðu0  u; X0  XÞ ¼ f ðn; gÞ

ð2Þ

where n ¼ X0  X is the relative position in the reference configuration, g ¼ u0  u denotes the relative displacement of two interacting material points. Moreover, the pairwise force function f is required to have the properties, i.e. f ðn; gÞ ¼ f ðn; gÞ, which assures conservation of linear momentum, and ðn þ gÞ  f ðn; gÞ ¼ 0, which assures conservation of angular momentum. For a micro-elastic material defined by Silling [21], the pairwise force function can be derived from scalar-valued micropotential function wðn; gÞ as

f ðn; gÞ ¼

@wðn; gÞ @g

ð3Þ

where wðn; gÞ is a micro-elastic strain energy, i.e. strain energy density per unit volume stored in the material. The micropotential function for the linear micro-elastic materials can be defined as

wðn; gÞ ¼

cðnÞs2 ðn; gÞjnj 2

ð4Þ

where cðnÞ designates the micro-modulus constant and is a constitutive parameter in the BB-PD model, and sðn; gÞ denotes the stretch of the bond which can be expressed as follows:

s¼

jn þ gj  jnj jnj

ð5Þ

where jnj is length of the bond in the reference configuration, and jn þ gj is the length of deformed bond in the deformed configuration. Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

5

Substituting Eq. (4) into Eq. (3), the pairwise force function through a PD bond for a homogeneous linear isotropic material is written as

f ¼ f ðn; gÞ ¼ cs

nþg jn þ gj

ð6Þ

3. Conjugated-bond linear elastic model 3.1. Motivation In the standard BB-PD model, only pairwise micro-potential function contains a micro-modulus constant, which is related to the function of the distance between the two interacting points. Thus, there only exists one mechanical parameter in the standard BB-PD model. However, in the macro scale, two macro mechanical parameters (Young’s modulus and Poisson’s ratio) govern the elastic systems. Therefore, the discrete elastic system described by pairwise micro-potential function in the standard BB-PD is restricted by a fixed Poisson’s ratio. In order to overcome the limitation of a fixed Poisson’s ratio in standard BB-PD, the micro-potential function, which is applied to describe the interaction between discrete particles in one influencing horizon, should be developed. The Stillinger-Weber (SW) potential functions [43,45] indicates that the bond energy is not only related to the bond stretch, but also related to the bond angle subtended by the given bond and other bonds [44,46]. Based on the Stillinger-Weber (SW) potential functions [43,45], the correlated lattice bond cell was proposed to simulated fracture problems by Zhang and Chen [44]. In the present study, a novel conjugated bond linear elastic constitutive model is developed in the framework of BB-PD based on the Stillinger-Weber (SW) potential functions [43–45]. Two micro-mechanical constants in the proposed conjugated BB-PD correspond to the two macro-mechanical constants in the conventional elastic mechanics, which overcome the limitation of the fixed Poisson’s ratio in the standard BB-PD. 3.2. Characterization of the strain energy In the novel conjugated bond model, interaction between two material points is described by the two- and three-point interaction. The energy of the bond nij is not only related to the deformation of its length, but also related to the bond angles subtended by itself and its conjugate bonds nik . For a horizon of a given particle Xi with N material points, the total number of bonds nij is N  1 and the total number of conjugated bonds nik is N  2, as shown in Fig. 2. Therefore, the strain energy density at a given material point in one horizon in the novel conjugated bond model can be expressed as PD W PD ¼ W PD l þ Wb

ð7Þ

where W PD and W PD l b are the stretch strain energy density and rotation strain energy density, respectively.

Fig. 2. Conjugated bond vectors subtended at a given material point Xi with N material points in the horizon of Xi .

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

6

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

The stretch and rotation strain energy density at a material point can be respectively calculated using the strain energy density induced by the two- and three-point interaction as

W PD l ¼

1 2

W PD b ¼

1 2

Z

Hd

wl ðnij ÞdV nij

Z

ð8Þ

Z

Hd

Hdh

wb ðnij ; nik ; hjik ÞdV nik dV nij

ð9Þ

jik

in which, wl ðnij Þ denotes the normal stretch strain energy micro-potential, which is related to the bond nij between material points Xi and Xj , wb ðnij ; nik ; hjik Þ designates the rotation strain energy micro-potential, which is related to the rotation between bond nij and bond nik . For the stretch of bond nij , length of bond between material points Xi and Xj in the reference configuration is jnij j ¼ Xj  Xi . However, in the deformed configuration, length of deformed bond can be expressed as jnij þ gij j ¼ ðXi þ ui Þ  ðXj  uj Þ. In the small deformation, according to Cauchy-Born rule [47], the lengths of the normalized bond nij can be determined as

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jnij þ gij j ¼ jnij j vT FT Fv ¼ jnij j vT ð2E þ IÞv

ð10Þ

in which, F is the deformation gradient, E is the Lagrangian strain tensor in the bond, I is the identity tensor,

v ¼ ½cos u; sin uT is the orientation vector of the bond in the reference configuration for 2D plane strain problems. Moreover, in the small deformation cases, the Lagrangian strain tensor is reduced to the Cauchy strain, i.e. E ! e. Then, the deformed bond length between material points can be obtained as

jnij þ gij j ¼ jnij jvm emn vn þ jnij j

ð11Þ

From Eq. (10), the strain energy micro-potential induced by the normalized stretch of bond nij in the linear elastic model can be expressed as:

wl ðnij Þ ¼

1 2 cn jnij jðvm emn vn Þ 2

ð12Þ

where cn is the normal stiffness of bonds, eij denotes the strain tensor in the bond nij .The macro-elastic stretch strain energy density due to the normal stretch of bond nij stored at material point Xi is obtained by integrating wl ðnij Þ over the horizon HX of Xi , i.e. integrating all bonds connected to Xi

W PD l ¼

1 2

Z

HX

wl ðnij ÞdV Xj

ð13Þ

i

Substituting Eq. (12) into Eq. (13), the PD macro-elastic energy density for a homogeneous linear elastic material is obtained as

W PD l ¼

1 2

Z 0

d

Z 2p 1 2 cn ðvm emn vn Þ jnj2 dudjnj 2 0

ð14Þ

where ðjnj; uÞ represents the integration in polar coordinates for 2D plane strain problems. For the rotation of a pair of conjugated bonds nij and nik , lengths of the conjugated bonds in the reference configuration can be expressed as

jnij j ¼ jXj  Xi j

ð15aÞ

jnik j ¼ jXk  Xi j

ð15bÞ

where Xi , Xj and Xk are coordinates of the three material points in the references configuration, as shown in Fig. 3. However, in the deformed configuration, lengths of the conjugated bonds can be respectively expressed as follows:

jnij þ gij j ¼ xj  xi ¼ jnij jvm emn vn þ jnij j

ð16aÞ

jnik þ gik j ¼ xk  xi ¼ jnik jwm emn wn þ jnik j

ð16bÞ

where xi , xj and xk are the coordinates of three material points in the deformed configuration, respectively. In the reference  ; sin u   are the unit orientation vectors of conjugated bonds nij and nik in the configuration, v ¼ ½cos u; sin u and w ¼ ½cos u polar coordinates, respectively. As shown in Fig. 3, the bond angle hjik between a pair of conjugated bonds nij and nik in the reference configuration and deformed configuration can be respectively written in the following form [44,46]: T

T

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

7

Fig. 3. Deformation of a pair of conjugated bonds.

0

1

0

1

vT w

vT w

B C B C hjik0 ¼ arccos @pffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiA ¼ arccos @pffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiA vT v  wT w vT v  wT w 0 hjik

ð17Þ

0

1

1

vT FT Fw vT ð2E þ IÞw B C B C ¼ arccos @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA ¼ arccos @pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA T T T T T T v ð2E þ IÞv  w ð2E þ IÞw v F Fv  w F Fw

ð18Þ

where F is the deformation gradient, E is the Lagrangian strain tensor, I is the identity tensor, v and w are the unit orientation vectors of conjugated bonds in the reference configuration. From Eqs. (17) and (18), the strain energy micro-potential induced by the rotation of a pair of conjugated bonds nij and nik in linear elastic model can be written as

wb ðnij ; nik ; hjik Þ ¼

1 kt ðhjik  hjik0 Þ2 2

ð19Þ

where kt designates the rotation resistance stiffness due to rotation of conjugated bonds, hjik0 denotes bond angles between nij and nik in the reference configuration, and hjik is bond angles between nij and nik in deformed configuration. From Eqs. (17) and (18), the bond angle is a continuous function of strain tensors, whose derivatives with respect to strain are extremely complicated. Bond angle of the conjugated bonds can be simplified by expanding bond angle into Taylor series:

hjik ðeÞ  hjik ð0Þ þ

@hjik ð0Þ @ 2 hjik ð0Þ 2 eþ  e þ  @e @ e2

ð20Þ

For simplicity, in this paper, the high strain gradient is neglected and will be studied in the future. Therefore, only the first-order term are considered to approximate the bond angle variation of the conjugated bonds, which is written as follows:

hjik ðeÞ  hjik ð0Þ þ

@hjik ð0Þ e @e

ð21Þ

where @hjik ð0Þ=@ e is the rate of rotation of conjugated bonds, which is labeled as Hjik in the following sections. Thus, the rotation angle of conjugated bonds in the deformed configuration is written as

Hjik  e ¼

@hjik ð0Þ  e  hjik ðeÞ  hjik ð0Þ @e

ð22Þ

Substituting Eqs. (17) and (18) into Eq. (22), the rotation angle of conjugated bonds in the deformed configuration can be re-written in the terms of strain tensors [44].

2vT ew þ vT w

A A11 

12 Hjik  e ¼ Q½e ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi

2v e v þ v v  T

T

2w ew þ w w T

T

A22

ð23Þ

where the scalars A11 , A22 and A12 are respectively written in the following forms.

A11 ¼ 2vm emn vn þ vm vm

ð24aÞ

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

8

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

A22 ¼ 2wm emn wn þ wm wm

ð24bÞ

A12 ¼ 2vm emn wn þ vm wm

ð24cÞ

Then, the first derivation of hjik ðeÞ with respect to

e is

  @hjik ðeÞ 1 @Q 1 1 A12 m n A12 m n ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2vm wn  v v  w w @ eij A11 A22 @e A11 A22 1  Q2 1  Q2

ð25Þ

Under the deformed state, according to Eqs. (21)–(24), the rate of rotation of conjugated bonds Hjik is written as :

Hij ¼

@hjik ð0Þ 1 ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ð2vm wn  vT w  vm vn  vT w  wm wn Þ @e 2 1  ðvT wÞ

ð26Þ

From Eqs. (22), (26) and (19), the strain energy micro-potential due to rotation of conjugated bonds can be determined by: 2

wb ðnij ; nik ; hijk Þ ¼ 12 kt ðHmn emn Þ  2 m n mn 1 m n T m n T ¼ 12 kt  pffiffiffiffiffiffiffiffiffiffiffiffiffiffi v w  v w  v v  v w  w w Þ e  ð2 2 T

ð27Þ

1ðv wÞ

For 2D plane strain problems, strain energy density in the novel conjugated bond model of BB-PD, which is induced by the rotation of bonds nij and nik , is expressed in the following form:

W PD b

1 ¼ 2

Z

d

Z 2p (Z 0

0

d

0

) Z 2p 1 1 m n mn 2 m n T m n T  kt  ½ð2v w  v w  v v  v w  w w Þe  jnjdudjnj du 2 1  ðvT wÞ2 0

ð28Þ

 are the angles of conjugated bonds in where jnj and j nj respectively represent the lengths of two conjugated bonds, u and u the polar coordinate. Therefore, for 2D plane strain problem, the two parts of strain energy density in the novel conjugated BB-PD can be written in the following form:

W PD l ¼

 cn pd3 h  2 3I1  4I2 48

ð29Þ

W PD b ¼

 kt p2 d2 h  2 I1  4I2 8

ð30Þ

where h is the thickness, I1 ¼ e11 þ e22 and I2 ¼ e11 e22  e212 , which implies that the strain energy density is objective. 3.3. Determination of conjugated bond parameters According to the characteristics of conjugated bond linear elastic model, the strain energy density around the horizon at a given material point can be written as





3 2 2 ¼ pd16hcn þ kt p8d h

PD W PD ¼ W PD l þ Wb    3   3 2 2 2 2 e211 þ e222 þ cn p24d h  kt p4d h e11 e22 þ pd12hcn þ kt p2d h e212

ð31Þ

For plan stain state, strain energy density obtained from classical elastic mechanics can be written as follows:

WC ¼

      1 2 2 K þ G e211 þ e222 þ K  G e11 e22 þ 2Ge212 2 3 3

ð32Þ

E E in which, K ¼ 3ð12 v Þ denotes the bulk elastic modulus, G ¼ 2ð1þv Þ designates shear modulus, E is Young’s modulus, and

v

is

Poisson’s ratio. By comparing strain energy density in the conjugated BB-PD model with that in the classical strain energy density for 2D plane strain state, the relationships between micro-parameters in the conjugated BB-PD model and macro-parameters in classical elastic mechanics can be expressed as follows:

pcn d3 h kt p2 d2 h 16

þ

8

¼

1 2 Kþ G 2 3

cn pd3 h kt p2 d2 h 2  ¼K G 24 4 3

ð33aÞ ð33bÞ

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

9

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

pd3 hcn 12

þ

k t p2 d 2 h ¼ 2G 2

ð33cÞ

From Eqs. (33a)–(33c), the relationships between the novel conjugated BB-PD parameters cn and kt can be written as:

cn ¼ kt ¼

6E ð1 þ v Þð1  2v Þpd3 h 2Eð1  4v Þ

ð1 þ v Þð1  2v Þp2 d2 h

ð34Þ ð35Þ

Based on Eqs. (34) and (35), two parameters are used in isotropy constitutive modeling to deduce the deformation under loads, whereas there is only one parameter, i.e. normal bond stiffness, in the standard BB-PD model. Some researchers developed some additional shear stiffness [38,48] and micro-polar model [49] in the BB-PD model. However, similar limitation of the Poisson’s ratio occurs in other BB-PD models. Moreover, this limitation of Poisson’s ratio is also inherent to the addition of a shear stiffness in lattice spring model [44,46,48,50]. However, in the present paper, there are two kinds of bond stiffness: normal bond stiffness and rotational resistant stiffness in the framework of the novel conjugated bond-based peridynamics, which overcomes the limitation of the fixed Poisson’s ratio in the standard bond-based peridynamics. In the present study, the rotational resistant stiffness kt is a generalized stiffness and can satisfy kt > 0, kt < 0 and kt ¼ 0, which is similar to the lattice model by Zhang and Ge [47]. It can be obtained from Eq. (34) that when the rotational resistant stiffness kt ¼ 0, the corresponding Poisson’s ratio is 0.25, which is equivalent to BB-PD. For the rotational resistant stiffness kt > 0, the corresponding Poisson’s ratio is lower than 0.25, which indicates that the tangential bond forces along the vertical direction of deformed bond have a strengthened effect on the normal bonds. In the case of rotational resistant stiffness kt < 0, the corresponding Poisson’s ratio is lower than 0.25, which indicates that the tangential bond forces have a weakened effect on the normal bonds. Therefore, the normal bond force density exerted on the point Xi in the deformed configuration can be written as follows:

f n ¼ c n  sn 

nþg jn þ gj  jnj  ^en ¼ cn  sn  ^en ¼ cn  jnj jn þ gj

ð36Þ

where cn is the normal stiffness of bonds, sn designates the normal stretch of the bond, ^ en is the unit orientation vector of the deformed bond nij , as shown in Fig. 4(a). The moment vector density exerted at the material point Xi , which is used to resist the variation of bond angle hjik subtended by the bonds nij and nik in Fig. 4(b), can be written in the following form:

mjik ¼ kt  ðhjik  hjik0 Þ

ð37Þ

where kt denotes the rotation resistance stiffness. Therefore, the tangent bond force density caused by rotation of a pair of conjugated bonds can be obtained as:

ft ¼

mjik kt  ðhjik  hjik0 Þ  ^et ¼  ^et jnij þ gij j jnij þ gij j

ð38Þ

where ^ et denotes the unit orientation vector which is along vertical direction of deformed bond nij , as shown in Fig. 4(c). The proof of the balance of angular moment is given in Appendix.

Fig. 4. Schematics of normal bond force, bond moment and tangent bond force.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

10

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Algorithm 1 is given to compute the internal forces at the material points with undamaged connecting bonds. Due to a loop over all the other bonds for a give material points, the computational cost for the novel conjugated bond-based peridynamics is larger than that for the standard bond-based peridynamics. Algorithm 1 Updating the internal forces 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

Do i = 1, totnode % Looping all the PD particles in the numerical model at one time step; f n ðiÞ 0 % Initiating the normal bond force density exerted on particle i; 0 % Initiating the tangent bond force density exerted on particle i; f t ðiÞ Do j = 1,numfam(i) % Looping the PD particles in the horizon of particle i; Labeling the number of particles in horizon i; Computing the distance between particle i and particle j in the reference configuration nij ; Computing the distance between particle i and particle j in the deformed configuration nij þ gij ;

Computing the normal force density f n ði; jÞ along the deformed bond nij þ gij based on Eq. (36); Summing the normal force density to compute the normal bond force density exerted on particle i, which is expressed as f n ðiÞ ¼ f n ðiÞ þ f n ði; jÞ; 0 % Initiating the tangent bond force density during the interaction between PD particle i and j; f t ði; jÞ Do k = 1,numfam(i) % Looping the PD particles in the horizon of particle i; Labeling the number of particles in horizon i; If (the number of particle k is equal to the number of particle j) then Stopping this loop; Elseif (the number of particle k is not equal to the number of particle j then Computing the distance between particle i and particle j in the reference configuration nik ; Computing the distance between particle i and particle j in the deformed configuration nik þ gik ; Computing the conjugated angle hjik0 in the reference configuration; Computing the conjugated angle hjik in the deformed configuration; Computing the variation of conjugated angle during deformation Dhjik ¼ hjik  hjik0 ; Summing over the N i  2 variations of conjugated angles; Enddo % Ending the loop of the particle k; Computing the moment vector density exerted at the particle i based on Eq. (37); Computing the tangent bond force density f t ði; jÞ based on Eq. (38); Summing the tangent bond force density exerted at the particle i f t ðiÞ ¼ f t ðiÞ þ f t ði; jÞ; Endodo % Ending the loop of the particle j; Enddo % Ending the loop of the particle j.

3.4. Damage model The method of damage initiation, propagation and coalescence in the proposed numerical model is through the permanent rupture of bonds. In the present study, an energy-based damage model is used to determine the rupture of bonds and simulate brittle fracture. If the strain energy density wn containing within a bond exceeds a certain critical value wc , the bond is assumed to be broken and the rupture of bonds is irreversible. Based on the theory of the conjugated BB-PD, the micropotential energy function in a bond in the conjugated BB-PD model can be expressed as: i X 1 2 1 cn s jnj þ kt ðDhjik Þ2 2 n 2 k¼1

N 2

wn ¼

ð39Þ

where cn and kt are respectively the normal stiffness and rotation resistant stiffness, sn denotes the normal stretch of a bond, Dhjik designates the variation of bond angle subtended by a pair of conjugated bonds,N i denotes the number of particles in the horizon of Xi . The critical strain energy density can be referred to the previous literature [23,51]. As shown in Fig. 5, when a crack passes through the horizon region of point A along the dash line ð0 < z < dÞ, every bond AB across the fracture surface has to be broken. If Gc is the fracture surface energy, i.e. the critical energy per unit fracture area required to create a crack surface, then the critical PD strain energy can be related to Gc . The critical strain energy density wc is associated with the bond connecting point A with the other point B. The sum of these energy density through integration is related to the fracture surface energy, which can be written as follows: Z d Z 2p Z d Z arccosðz=nÞ Gc ¼ wc n2 sin hdhdndudz ð40Þ 0

0

z

0

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

11

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 5. Evaluation of critical strain energy density for a material point.

By solving the integration, we can obtain the critical strain energy potential in the bond between two interacting material points as

wc ¼

4Gc

ð41Þ

pd4

When the above damage model is numerically implemented into the explicit dynamic code,wn in every bond at each time step is compared with the critical strain energy potential. If wn exceeds wc , the breakage of the bond n occurs. A historydependent scalar-valued function lðwn ; tÞ is used to describe the breakage of bonds. The failure criterion is expressed as follows:



lðwn ; tÞ ¼

0

wn > wc

1 otherwise

ð42Þ

where wc is the critical strain energy potential, which can be obtained from Eq. (41).For damage definition convenience, the damage at a point X at time t could be defined as:

R

DðX; tÞ ¼ 1 

HX

lðwn ; tÞdV R

HX

dV

ð43Þ

where DðX; tÞ is the damage value at a point X at time t, which varies from 0 to 1. When the damage value at point X is equal to one, it indicates the location of cracks. However, when the damage value at point X is smaller than one, it means the damage around cracks. The detail process of determining bond rupture criteria is presented in Algorithm 2. Algorithm 2 Cracking tracking algorithm 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.

Do i = 1, totnode % Looping all the PD particles in the numerical model at one time step; Do j = 1,numfam(i) % Looping the PD particles in the horizon of particle i; Computing the normal stretch bond strain energy potential based on Eq. (12); Computing the rotation strain energy potential based on Eq. (19); Computing the micro-potential wnij between PD particles i and j based on Eq. (39) If (micro-potential wnij > critical strain energy potential wc ) then Breakage of the bond nij occurs based on Eq. (42); Else The bond nij is intact based on Eq. (42); Endif Enddo Computing the damage of a given material point I based on Eq. (43); Enddo

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

12

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

3.5. Time stepping In the reference configuration, a body is uniformly discretized into particles with a certain volume [23]. All particles form a grid allowing deformation without restriction, which is different from FEM. After discretization, the integral-differential equation of motion in the conjugated BB-PD model can be formulated as

€ ni u

q

¼

N i 1 X

f n DV j þ

j¼1

N i 1 X

N i 2 X

j¼1

k1

! n

f t DV j þ bi

ð44Þ

€ ni denotes the acceleration on the particle Xi at time step n, f n designates the normal force density due to stretch of in which u n bond nij , f t represents the tangent force density due to the rotation of a pair of conjugated bonds nij and nik , and bi is the body force density exerted on the particle Xi at time step n. The displacement of particle Xi at time step ðn þ 1Þ can be obtained by approximating the acceleration in Eq. (43) in an explicit formation

€ ni ¼ u

unþ1 i

unþ1  2uni þ un1 i i Dt 2 ¼

"N 1 i Dt 2 X

q

f n DV j þ

j¼1

ð45Þ N i 1 X

N i 2 X

j¼1

k1

! f t DV j þ

# n bi

þ 2uni  un1 i

ð46Þ

where Dt designates the time step size, f n and f t respectively denote the normal bond forces and tangential bond forces n exerted on the particle Xi by other particles Xj within the horizon of the given particle Xi , bi is the body force density on the th particle Xi at the n time step. According to works by Silling and Askari [23], the stable time step should be satisfied by the following condition

Dt <

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2q PN i j¼1 C ij V j

ð47Þ

where C ij ¼ jCðX j  X i Þj ¼ j@f =@ gj is a scalar-valued function for micromodulus, V j denotes the volume of the particle Xj in the horizon of a given particle Xi .

Fig. 6. Geometry of a plant with a circular cutout under tensile loads.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

13

Fig. 7. Horizontal and vertical displacements of the numerical sample obtained from the proposed numerical method (a)–(b) and FEM (c)–(d).

4. A benchmark example 4.1. Numerical model In this section, an isotropic plate with a circular cutout subjected to a slow rate loads along the horizontal edges, which represents the quasi-static loads. The geometric and loading conditions are shown in Fig. 6(a), and the numerical sample is plotted in Fig. 6(b). The material properties are listed as follows: Young’s modulus is E ¼ 192 GPa, Poisson’s ratio is v ¼ 0:2, mass density is q ¼ 8000 kg=m3 . The tensile loading rate is u_ y ¼ 2:7541  107 m=s. In the numerical sample, 100  100 ¼ 10; 000 PD particles represent the plate with a circular cutout, and the time step is adopted as Dt ¼ 1:0  107 s.

4.2. Numerical results Fig. 7(a) and (b) shows the horizontal displacement field and vertical displacement field at the real time of 1000 s in the absence of failure. Compared with the numerical results obtained from FEM in Fig. 7(c) and (d), it can be found that the predicted horizontal and vertical displacement fields obtained from the proposed numerical method are in good agreement with those obtained from FEM. In addition, the variations of horizontal and vertical displacements along x-axis and y-axis obtained from the conjugated BB-PD and FEM are respectively plotted in Fig. 8. A close agreement can be observed among the conjugated BB-PD predictions and FEM results. It indicates the accuracy of numerical results obtained from the conjugated BB-PD codes. Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

14

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 8. Comparison of the horizontal displacement (a) and vertical displacement (b) along the central axes obtained from the conjugated BB-PD and the finite element method (FEM).

5. Numerical examples and results 5.1. Notched beam subjected to impact loads In order to examine the application of the proposed numerical method to simulate dynamic fracture growth, a concrete beam with an offset pre-existing crack subjected to the impact loads was investigated using the experimental method [52,53] and other numerical method [44,54,55]. The dimension of the edge notched beam is plotted in Fig. 9. According to works by Zhang and Chen [44] and Belytschko and Tabbara [54], the beam is subjected to the impact velocity at the upper edge in Fig. 9, which can be expressed as follows.

VðtÞ ¼

V 0 t=t 1

for t 6 t1

V0

for t > t 1

ð48Þ

in which t 1 ¼ 1:96  104 s and V 0 ¼ 0:065 m=s. In order to investigate the influence of mechanical parameters of materials on dynamic fracture propagation, three cases of edge notched beams with different mechanical parameters are simulated, as shown in Table 1. The other mechanical parameters are listed as follows: Poisson’s ratio is v ¼ 0:2 and material density is q ¼ 2400 kg=m3 . pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi From the shear wave velocity cs ¼ E=ð2ð1 þ v ÞqÞ, the Rayleigh wave velocity in case I can be estimated as cR  0:9cs . TherePlease cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

15

Fig. 9. An edge notched beam subjected to impact loads (a) geometric and loading conditions; (b) numerical discrete particles scheme.

Table 1 Different mechanical parameters of materials for three edge notched beams. Young’s modulus Case I

E ¼ 31:37 GPa

Case II

E ¼ 50:0 GPa

Case III

E ¼ 31:37 GPa

Critical stress intensity factor K cI ¼ 0:8 MPa m1=2 K cI ¼ 0:8MPa m1=2 K cI ¼ 1:6 MPa m1=2

fore, the Rayleigh wave velocity in the present numerical simulation is cR ¼ 2100 m=s. In the numerical sample, 228  76 ¼ 17; 328 PD particles represent the numerical notched beams. The spacing of PD particles with uniform distribution is taken as Dx ¼ 1 mm and the radius of the influence horizon is adopted as d ¼ 3Dx. As shown in Fig. 10, the fields of strain energy density during the crack growth process in Case I are plotted in the left column, and crack damage maps during the crack growth and deformation process in Case I are presented in the right column. It can be observed from Fig. 10(a) that the initiation of the crack occurs at 340 ls due to the concentration of strain energy density around the tip of the pre-existing notch. As the impact load increases, the crack continues to propagate towards the loading point in Fig. 10(b)–(e). Meanwhile, the strain energy density is concentrated around the propagating crack tip from 376 ls to 492 ls. Finally, the notched beam subjected to impact loads is completely failed at 540 ls. It can be observed from Fig. 11(a)–(c) that the final predicted crack trajectories in numerical samples in Case I, Case II and Case III are similar. Compared the final predicted crack trajectories obtained from the proposed numerical method with those obtained from the previous experimental observations [52] and the previous numerical results [44] in Fig. 11 (d) and (e), the present numerical results are in good agreement with the previous experimental observations and numerical results. Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

16

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 10. Fracture propagation process of elastic strain energy density (left column) and damage maps (crack trajectories, where the displacements of particles are magnified by a factor of 100 in the right column) in Case I at different times.

In addition, crack propagation speeds in Case I, Case II and Case III are plotted in Fig. 12. The profiles of crack propagation speed in the three cases are similar. The initiation of the crack in Case I occurs at 340 ls, and the initiation of the crack in Case II appears at 180 ls. It implies that the bigger Young’s modulus is, the more earlier the crack initiates. The maximum crack propagation speed in Case I is 855 m=s, which is equal to 0:41cR . While the maximum crack propagation speed in Case II is 727 m=s, which is equal to 0:346cR . By comparing the maximum crack propagation speeds in Case I and Case II, it is

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

17

Fig. 11. Comparison of the final predicted crack trajectories in the edge notched beams subjected to impact loads in (a) Case I, (b) Case II, (c) Case III obtained from the proposed numerical method, (d) the previous experimental observation [50] and (e) the previous numerical results [42].

Fig. 12. Crack growth speed versus time in cases I, II and III.

found that the bigger Young’s modulus is, the lower the maximum crack propagation speed is. Similar conclusions were drawn using corrected lattice bond cell by Zhang and Chen [44].

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

18

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 13. Effect of Poisson’s ratio on final predicted crack trajectories, where the displacements of PD particles are magnified by 100 times: (a) v ¼ 0:15, (b) v ¼ 0:25 and (c) v ¼ 0:35.

Fig. 14. Effect of Poisson’s ratio on the curves of crack speed versus time.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

19

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx Table 2 PD particle spacings and sizes of influence horizon. Label

PD particle spacing Dx (mm)

Size of the influence horizon d (mm)

Number of PD particles in the numerical sample

A B C

2 1 0.5

6 3 1.5

4332 17,328 69,312

By comparing crack propagation speeds in Case I and Case III in Fig. 12, it can also be observed that the larger the critical stress intensity factor is, the later the crack initiates, and the lower the maximum crack propagation speed is. In order to study the effects of Poisson’s ratio on crack propagation in the notched beam under impact loads, three numerical samples with Poisson’s ratio of 0:15, 0:25 and 0:35 in Case I are simulated, respectively. The final predicted crack propagation patterns in three numerical samples are plotted in Fig. 13. It can be found from Fig. 13 that Poisson’s ratio obviously affects the roughness of crack surface. The roughness of crack surface becomes more coarse with increasing Poisson’s ratio, as shown in Fig. 13(a)–(c). In addition, crack propagation speeds in the above three numerical notched beams are plotted in Fig. 14. It can be observed from Fig. 14 that cracks in the above three numerical samples with different Poisson’s ratios initiate at the almost same time. The maximum value of crack propagation speed becomes larger as Poisson’s ratio increases. For the numerical sample with Poisson’s ratio of 0:15, the maximum value of crack propagation speed is equal to 964 m=s, which is about 0:46cR . For the numerical sample with Poisson’s ratio of 0:25, the maximum value of crack propagation speed is about 804 m=s, which is equal to 0:38cR . For the numerical sample with Poisson’s ratio of 0:35, the maximum crack propagation speed is approximate 705 m=s, which equals 0:34cR . It can be reasonable to deduce that the declining tendency of the maximum crack propagation speed may be caused by the larger roughness of crack surfaces which is induced by the increase of Poisson’s ratio. According to works by Bobaru et al. [29], the numerical convergence of PD results can be divided into two types: the mconvergence and the d-convergence. For the m-convergence, the size of influence horizon d is fixed and the number of particles in one horizon m increases. The numerical PD approximation converges to the exact non-local PD solution for the given d. While, for the d-convergence, the number of particles in one horizon m is fixed and the size of influence horizon d decreases. In the d-convergence, the numerical PD result converges to an approximation of the classical solution. Moreover, the larger m is, the closer this approximation becomes. In this present numerical example, in order to investigate the dconvergence of the conjugated BB-PD, three numerical edge notched beams subjected to impact loads are simulated. In addition, the m-convergence of the conjugated BB-PD will be discussed in the next section. For the fixed m of 3, the effect of the various d on the numerical results is investigated by using three numerical edge notched beams respectively marked as A, B and C. The geometric dimensions and mechanical parameters in the present numerical samples are the same as those in Case I. Various PD particle spacings and sizes of influence horizon in the three numerical beams are listed in Table 2. The final predicted crack propagation patterns in the above three numerical samples are plotted in Fig. 15(a)–(c). It can be found from Fig. 15(a)–(c) that the predicted failure patterns of the notched beams subjected to impact loads are similar, which indicates the numerical convergence of the proposed numerical method. However, the width of the predicted failure zone in the numerical sample increases with increasing the influence horizon d. The profiles of crack propagation speeds in three numerical samples are plotted in Fig. 16(d). It can be obviously found from Fig. 16(d) that the crack propagation speeds in three numerical samples are similar and the maximum crack propagation speed slightly increases with an increase in horizon size d.

5.2. Crack branching Crack branching problem is a typical dynamic fracture problem, which was studied by some researchers to verify different numerical methods [1,11,28,32,34,44,57,58]. The schematic of the experiment is plotted in Fig. 16. Labortory experiments on specimens with similar dimensions were conducted by Ramulu and Kobayashi [59], Ravi-Chandar [60], Sharon et al. [61], Sharon and Fineberg [62] and Fineberg et al. [63]. The material properties are listed as follows: q ¼ 2450 kg=m3 , E ¼ 32 GPa and v ¼ 0:20. From works by Song et al. [64], the energy release rate is Gc ¼ 3 N=m. The dilatational wave speed is cd ¼ 3809:5 m=s and the shear wave speed is cs ¼ 2332:8 m=s. The corresponding Rayleigh wave speed is cR ¼ 2119:0 m=s. It can be seen from Fig. 16 that a tensile stress ry ¼ 1 MPa is applied on the top and bottom boundaries and the time history of the load is a step function, which is expressed as follows:



ry ðtÞ ¼ in which

r0 t=t0 for t 6 t0 r0 for t > t 0

ð49Þ

r0 ¼ 1:0 MPa and t0 ¼ 50 ls.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

20

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 15. Final predicted crack trajectories in the numerical samples with the fixed m for the different influence horizon sizes at 540 ls: (a) notched beam A, (b) notched beam B, (c) notched beam C, and (d) the d-convergence for crack propagation speed versus time profiles.

In this numerical simulation, 17; 200 PD particles represent the rectangular plate. The spacing of PD particles is equal to Dx ¼ 0:0005 m and the radius of influence horizon is taken as d ¼ 0:0015 m. Moreover, the time step is taken as Dt ¼ 0:05 ls, which is smaller than the critical time step to obtain the stable solutions from Eq. (46). Fig. 17 shows the crack evolution predicted by the proposed numerical method at different time. It can be observed from Fig. 17 that crack initiation occurs at around 18 ls due to the concentration of strain energy density around the pre-existing crack tip. As the tensile loads increase, the typical mode I crack propagates along the horizontal direction with minor branches starting from the main crack, but only extending to a small length. At around 40 ls, the main crack is splitted into two main branches, which grow until they reach the right boundary of the numerical sample at 62 ls. By comparing with the previous experimental observations and other numerical results in Fig. 18, it can be found that the predicted crack branching pattern obtained from the proposed numerical method is in good agreement with those obtained from other experimental and numerical results. In addition, the profile of crack propagation speed versus time is plotted in Fig. 19. The crack propagation speed is in good agreement with other numerical results, such as CZM by Song et al. [64], DM by Braun and Fernández-Sáez [58], corrected lattice bond cell by Zhang and Chen [44], and etc. The initiation of the mode I crack at around 18 ls obtained from the conjugated BB-PD is in good agreement with that obtained from CZM, DM and corrected lattice bond cell. As shown in Fig. 19, crack branching occurs at around 38 ls. It is found from Fig. 19 that the crack propagation speed suddenly drops just before the occurrence of crack branching, which is consistent with the prediction by Zhang and Chen [44] and Zhou et al. [66]. Moreover, Fig. 19 shows that crack propagation speed fluctuates in an instable range, which does not exceed the Rayleigh wave speed. In order to investigate the m-convergence of the conjugated BB-PD, the numerical results of the crack-path and crack propagation speed for the dynamic crack branching problem [34] are discussed. Three groups of numerical samples, which are respectively labeled as ➀, ➁ and ➂, are simulated using the proposed numerical method. The size of the influence horizon, the spacing of material PD particles, the value of mð¼ d=DxÞ and the number of PD particles in the numerical sample are listed in Table 3. For the fixed horizon d of 1:5 mm, the m-convergence study is performed on the above three numerical samples. For the coarse numerical sample ➀ with the spacing Dx of 1:0 mm, there are 9 PD particles in the horizon ðm ¼ 1:5Þ. For the numerical sample ➁ with the spacing Dx of 0:5 mm, there are 25 PD particles in the horizon ðm ¼ 3Þ. For the numerical sample ➂ with the spacing Dx of 0:25 mm, there are 81 PD particles in the horizon ðm ¼ 6Þ.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

21

Fig. 15 (continued)

Fig. 16. A plate with an edge pre-existing crack.

Fig. 20(a)–(c) shows the final predicted crack branching patterns in the numerical samples with different values of m. Fig. 20(d) presents the crack propagation speed for three influence horizons. It can be observed from Fig. 20(a)–(c) that crack branching patterns in all cases are similar. As shown in Fig. 20(a)–(d), the initiation of cracks in all numerical samples occurs at around 18 ls, and crack branching appears at the location of approximate 0:069 m away from the left-side of the plate at

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

22

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 17. Fracturing process of crack branching: (a) strain energy density and (b) damage maps.

around 38 ls. It can also be found from Fig. 20(a)–(c) that the crack branching patterns become smoother and some small crack branches can be observed more easily with increasing m. Moreover, the curves of crack propagation speed versus time in three numerical samples in Fig. 20(d) are similar. Crack branching occurs after the crack propagations speed reaches

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

23

Fig. 18. Comparison of final crack trajectories obtained from (a) the proposed numerical method and (b) the previous experimental observation by Ramulu and Kobayashi [57] and (c)–(f) other numerical methods by Wildman and Gazonas [65], Ha and Bobaru [32], Song et al. [62], Braun and Fernández-Sáez [56].

Fig. 19. Comparison of crack growth velocity predicted by different numerical methods.

0:53cR  0:69cR , which is equal to 1123—1462 m=s. The crack propagation speeds in three numerical samples unstably fluctuate in an certain range, which do not exceed the Rayleigh wave speed. The crack branching patterns and crack propagation speeds with different horizons in Fig. 20 indicate the m-convergence. In order to investigate the effect of PD particle arrangement on crack branching, two numerical samples with uniform particle arrangement and hexagon particle arrangement are respectively conducted using the proposed numerical samples.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

24

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Table 3 Parameters for m-convergence in the conjugated BB-PD. Sample

Size of the influence horizon (mm)

Spacing of material PD particles (mm)

mð¼ d=DxÞ

The number of PD particles in the numerical sample

➀ ➁ ➂

1.5 1.5 1.5

1.0 0.5 0.25

1.5 3 6

4600 17,200 66,400

Fig. 20. Final predicted crack trajectories for the m-convergence: (a) Numerical sample ➀; (b) Numerical sample ➁; Numerical sample ➂, and (d) the curves of crack propagation speed versus time.

Fig. 21 shows crack branching and deforming processes of two numerical samples. It can be observed from Fig. 21 that crack initiation and branching in the numerical sample with hexagon particle arrangement occurs much earlier than those in the numerical sample with uniform particle arrangement. PD particle arrangement has slight influence on the time of crack initiation and branching. In addition, the ‘mirror-mist-hackle’ phenomenon can be observed from the crack branching and deforming processes of two numerical samples, which indicates that the dynamic crack propagation successively experiences smooth, rough and zigzag stage. The numerical reproduction of this typical phenomenon obtained from the proposed numerical method is similar to those obtained by Abraham et al. [67,68], Fineberg et al. [69], Cramer et al. [70], Ravi-Chandar and Knauss [60], Yang et al. [71] and Zhang and Chen [44].

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

25

Fig. 20 (continued)

5.3. Fragmentation of a cylinder under internal pressure According to works by Rabczuk and Belytschko [11], fragmentation of a cylinder subjected to internal pressure is modeled using the proposed numerical method. The inner radius of the cylinder is 80 mm and the outer radius is 150 mm. The mechanical parameters are listed as follows: Young’s modulus E ¼ 200 GPa, the density q ¼ 7800 kg=m3 , Poisson ratio v ¼ 0:2 and the critical stress intensity factor K cI ¼ 6:0 GPa m1=2 . A pressure of p ¼ p0 et=t0 is applied with p0 ¼ 10; 000 MPa and t 0 ¼ 0:1 ms (t is the real time). The numerical sample containing 12,668 discrete material PD particles, which are uniformly distributed with orientation of 0 , are simulated, as shown in Fig. 22(a). Moreover, in order to investigate the effects of discrete PD particle orientation on the predicted crack paths, other three numerical samples with discrete PD particle orientation of 15 , 30 and 45 in Fig. 22(b)–(d) are respectively conducted using the proposed numerical method. Figs. 23–26 show the fracturing and deformation processes of the numerical cylinders with different PD particle orientations under internal pressure. It can be observed from Fig. 23(a) that the tensile macro-cracks alternately occur around the inner wall of the cylinder subjected to internal pressure at 50 ls. When internal pressure increases, the tensile macro-cracks continue to propagate along the radial directions. Finally, the cylinder breaks into several fragments and some damaged PD particles fly apart at 75 ls.The numerical results obtained from the proposed numerical method is in good agreement with those obtained by Rabczuk and Belytschko [11].

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

26

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 21. The effect of the discrete PD particle arrangements on crack branching patterns which are magnified by 100 times: (a) uniform distribution and (b) hexagonal distribution.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

27

Fig. 22. Numerical samples with the different PD particle orientations of (a) 0 ; (b) 15 ; (c) 30 and (d) 45 .

Fig. 23. Fracturing and deformation process of a cylinder under internal pressure at different times: (a) 50 ls; (b) 60 ls; (c) 70 ls and (d) 75 ls, respectively.

Similar fracturing processes and final predicted fragments can also be observed from Figs. 24–26 in the numerical cylinders with different PD particle orientations. It can be obtained from numerical results that PD particle orientation slightly affects fragmentation patterns. However, the fragmenting patterns of numerical cylinders with different PD orientations are not the same, i.e. the number of fragments. The number of fragments for different PD particle orientations is listed in Table 4. It can be found from Table 4 that the number of fragments decreases when the PD particle orientation varies from 0 to 45 .

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

28

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 24. Fracturing and deformation process of a cylinder with discrete PD particle orientation of 15 at different times: (a) 50 ls; (b) 60 ls; (c) 70 ls and (d) 75 ls.

6. Conclusions A novel conjugated bond linear elastic model is proposed in the framework of the standard bond-based peridynamics (BB-PD). Based on the modified SW potential, micro-potential function of the bond between two interacting points in the proposed numerical method consists of the normal stretch strain energy micro-potential and rotation strain energy micro-potential. In this paper, the micro-potential function of a bond is in the linear elastic case, which is related to the two- and three-point potential. The micro-mechanical parameters, i.e. the normal bond stiffness and rotation resistance stiffness, in the conjugated BB-PD are derived and calibrated using strain energy equivalence to a macroscopic classical continuum model, which overcomes limitation of the fixed Poisson ratio in the standard BB-PD. In addition, an energy-based damage model is implemented into the conjugated BB-PD model, which establishes the relationship between the proposed numerical method and linear elastic fracture mechanics. By comparing with the numerical results obtained from FEM, a benchmark example is performed to prove the ability and accuracy of the proposed numerical method. Three numerical dynamic fracture problems are simulated using the proposed numerical method. Some characteristics of the conjugated BB-PD are investigated and discussed. The conclusions are drawn as follows: (1) The present numerical results for the dynamic fracture problems are in good agreement with the previous experimental observations and numerical results, which indicates that the proposed numerical method is efficient for dynamic fracture problems. (2) Poisson’s ratio of materials not only affects the maximum crack propagation speed, but also influences the roughness of the new crack face. (3) d-convergence and m-convergence in the conjugated BB-PD are respectively validated by the notched beams subjected to impact loads and crack branching. (4) The effects of uniform particle arrangement and hexagon particle arrangement on crack branching are studied. The final predicted crack trajectories are similar but the roughness of the new crack face is slightly affected by the discrete PD particle arrangements.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

29

Fig. 25. Fracturing and deformation process of a cylinder with discrete PD particle orientation of 30 at different times: (a) 50 ls; (b) 60 ls; (c) 70 ls and (d) 75 ls.

Fig. 26. Fracturing and deformation process of a cylinder with discrete PD particle orientation of 45 at different times: (a) 50 ls; (b) 60 ls; (c) 70 ls and (d) 75 ls.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

30

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx Table 4 The dependence of the number of fragmentations of the cylinder under internal pressure on the orientation of PD particles. Orientation of PD particles (°)

The number of fragmentations

0 15 30 45

16 12 12 11

(5) The effects of discrete PD particle orientation in the numerical cylinder under internal pressure are investigated. Final fragmentations of the numerical samples with different discrete PD particle orientation are similar. When the PD particle orientation varies from 0 to 45 , the number of fragments decreases, which implies the slight effects of PD particle orientations on crack propagation patterns in the dynamic fracture problems. Due to the use of multiple bonds in the proposed numerical method, the computational cost of the proposed numerical method increases. Furthermore, the proposed numerical method will be extended to three dimension, and the GPU-based parallel programming technique will also be implemented into the proposed numerical method to solve fracturing problems. Acknowledgement The present work is supported by National Natural Science Foundation of China (Grant Nos.51325903 and 51679017), Graduate Scientific Research and Innovation foundation of Chongqing, China (Grant No. CYB16012), National Basic Research Program of China (973 Program: 2014CB046903), and Natural Science Foundation Project of CQ-CSTC (Grant Nos. cstc2015jcyjys30006, cstc2016jcyjys0005 and cstc2017jcyjA1250), which are gratefully acknowledgement. Appendix A The direct normal bond force induced by bond nij at point Xi within horizon HXi is written as

1 f n ðnij ; gij Þ  DV Xi  DV Xj 2

8Xj 2 HXi

ðA:1Þ

The reaction normal bond force induced by bond nji at point Xi within horizon HXj is

1 1  f n ðnji ; gji Þ  DV Xj  DV Xi ¼  f n ðnij ; gij Þ  DV Xj  DV Xi 2 2

8Xi 2 HXj

ðA:2Þ

Similarly, the direct tangent bond force induced by conjugated bond angle hjik acting at the point Xi within horizon HXi is obtained as

1 f t ðnij þ gij ; nik þ gik Þ  DV Xi  DV Xj 2

8Xj ; Xk 2 HXi

ðA:3Þ

The reaction tangent bond force induced by conjugated bond angle hijm exerted at the point Xi within horizon HXj is expressed as

1  f t ðnji þ gji ; njm þ gjm Þ  DV Xi  DV Xj 2

8Xi ; Xm 2 HXj

ðA:4Þ

The angular momentum H0 of a fixed set of material points at time t in the bounded body X is given by

Z

H0 ¼

X

_ i ; tÞdX ðXi þ ui Þ  qðXi ÞuðX

ðA:5Þ

While the total torque P0 about the origin is given as follows:

R P0 ¼ X ðXi þ ui Þ  bðXi ; tÞdX R R R R þ X HX ðXi þ ui Þ  12 f n ðnij ; gij ÞdV Xj dX  X HX ðXi þ ui Þ  12 f n ðnji ; gji ÞdV Xi dX i j R R R R þ X HX ðXi þ ui Þ  12 f n ðnik ; gik ÞdV Xk dX  X HX ðXi þ ui Þ  12 f n ðnki ; gki ÞdV Xi dX R R i Rk R þ X HX ðXi þ ui Þ  12 f t ðnij ; gij ; nik ; gik ÞdV Xj dX  X HX ðXi þ ui Þ  12 f t ðnji ; gji ; njm ; gjm ÞdV Xi dX i j R R R R þ X HX ðXi þ ui Þ  12 f t ðnik ; gik ; nij ; gij ÞdV Xk dX  X HX ðXi þ ui Þ  12 f t ðnki ; gki ; nkl ; gkl ÞdV Xi dX i

ðA:6Þ

k

_ 0 ¼ P0 results in Thus, the balance of angular momentum H

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

R

€ ðXi ; tÞ  bðXi ; tÞdX ¼ þ ui Þ  ½qðXi Þu R R R R þ X HX ðXi þ ui Þ  12 f n ðnij ; gij ÞdV Xj dX  X HX ðXi þ ui Þ  12 f n ðnji ; gji ÞdV Xi dX i j R R R R þ X HX ðXi þ ui Þ  12 f n ðnik ; gik ÞdV Xk dX  X HX ðXi þ ui Þ  12 f n ðnki ; gki ÞdV Xi dX R R i Rk R þ X HX ðXi þ ui Þ  12 f t ðnij ; gij ; nik ; gik ÞdV Xj dX  X HX ðXi þ ui Þ  12 f t ðnji ; gji ; njm ; gjm ÞdV Xi dX i j R R R R þ X HX ðXi þ ui Þ  12 f t ðnik ; gik ; nij ; gij ÞdV Xk dX  X HX ðXi þ ui Þ  12 f t ðnki ; gki ; nkl ; gkl ÞdV Xi dX

31

X ðXi

i

ðA:7Þ

k

Because f n ðnij ; gij Þ ¼ f n ðnji ; gji Þ ¼ 0 for Xj R HXi and f n ðnik ; gik Þ ¼ f n ðnki ; gki Þ ¼ 0 for Xk R HXi , terms on RHS of Eq. (A.7) can be rewritten to include all of the material points in the body X, which is given by

R

€ ðXi ; tÞ  bðXi ; tÞdX ¼ þ ui Þ  ½qðXi Þu R R R R 1 ðX þ u Þ  f ðn ; gij ÞdXdX  X X ðXi þ ui Þ  12 f n ðnji ; gji ÞdXdX i i X X 2 n ij R R R R þ X X ðXi þ ui Þ  12 f n ðnik ; gik ÞdXdX  X X ðXi þ ui Þ  12 f n ðnki ; gki ÞdXdX R R R R þ X HX ðXi þ ui Þ  12 f t ðnij ; gij ; nik ; gik ÞdV Xj dX  X HX ðXi þ ui Þ  12 f t ðnji ; gji ; njm ; gjm ÞdV Xi dX i j R R R R þ X HX ðXi þ ui Þ  12 f t ðnik ; gik ; nij ; gij ÞdV Xk dX  X HX ðXi þ ui Þ  12 f t ðnki ; gki ; nkl ; gkl ÞdV Xi dX

X ðXi

þ

i

ðA:8Þ

k

Then, by exchanging the dummy variables as Xi $ Xj and ui $ uj in the first term on RHS of Eq. (A.8), and Xi $ Xk and ui $ uk in the third term on RHS of Eq. (A.8), the second and fourth terms on RHS of Eq. (A.8) become

Z Z

Z Z

X

X

1 ðXi þ ui Þ  f n ðnij ; gij ÞdXdX ¼ 2 X

1 ðXj þ uj Þ  f n ðnji ; gji ÞdXdX 2 X

Z Z

Z Z

X

X

1 ðXi þ ui Þ  f n ðnik ; gik ÞdXdX ¼ 2 X

1 ðXk þ uk Þ  f n ðnki ; gki ÞdXdX 2 X

ðA:9aÞ ðA:9bÞ

Thus, Eq. (A.8) can be rewritten in the following forms

R

€ ðXi ; tÞ  bðXi ; tÞdX ¼ þ ui Þ  ½qðXi Þu R R ½ðX þ u Þ  ðXi þ ui Þ  12 f n ðnji ; gji ÞdXdX j j X X R R þ X X ½ðXk þ uk Þ  ðXi þ ui Þ  12 f n ðnki ; gki ÞdXdX R R R R þ X HX ðXi þ ui Þ  12 f t ðnij ; gij ; nik ; gik ÞdV Xj dX  X HX ðXi þ ui Þ  12 f t ðnji ; gji ; njm ; gjm ÞdV Xi dX i j R R R R þ X HX ðXi þ ui Þ  12 f t ðnik ; gik ; nij ; gij ÞdV Xk dX  X HX ðXi þ ui Þ  12 f t ðnki ; gki ; nkl ; gkl ÞdV Xi dX

X ðXi

þ

i

ðA:10Þ

k

in which the first and second terms on RHS of Eq. (A.10) are equal to zero, respectively. For the tangent bond force part, the third term on RHS of Eq. (A.10) can be recast by exchanging ðXi þ ui Þ  12 f t ðnij ; gij ; nik ; gik ÞXi $ Xj , Xk $ Xm , ui $ uj and uk $ um as follows:

Z Z X

HX

i

1 ðXi þ ui Þ  f t ðnij ; gij ; nik ; gik ÞdV Xj dX ¼ 2

Z Z X

HX

j

1 ðXj þ uj Þ  f t ðnji ; gji ; njm ; gjm ÞdV Xi dX 2

ðA:11aÞ

Similarly, by exchanging Xi $ Xk , Xj $ Xl , ui $ uk and uj $ ul as follows:

Z Z X

HX

i

1 ðXi þ ui Þ  f t ðnik ; gik ; nij ; gij ÞdV Xk dX ¼¼ 2

Z Z X

HX

j

1 ðXk þ uk Þ  f t ðnki ; gki ; nkl ; gkl ÞdV Xk dX 2

ðA:11bÞ

Thus, Eq. (A.10) can be re-written in the following forms:

R

€ ðXi ; tÞ  bðXi ; tÞdX ¼ þ ui Þ  ½qðXi Þu R R 1 ½ðX þ u Þ j j  ðXi þ ui Þ  2 f t ðnji ; gji ; njm ; gjm ÞdV Xi dX X HX j R R þ X HX ½ðXk þ uk Þ  ðXi þ ui Þ  12 f t ðnki ; gki ; nkl ; gkl ÞdV Xi dX R R k R R ¼ X HX ðnij þ gij Þ  12 f t ðnji ; gji ; njm ; gjm ÞdV Xi dX þ X HX ðnik þ gik Þ  12 f t ðnki ; gki ; nkl ; gkl ÞdV Xi dX

X ðXi

þ

j

ðA:12Þ

k

When we consider Eqs (37) and (38), it can be obtained that magnitude of moment vector induced by deformed bond ðnij þ gij Þ and tangent bond force density is the same as that caused by deformed bond ðnik þ gik Þ and tangent bond force density with opposite directions. Therefore, the balance of angular momentum can be satisfied as

Z

X

€ ðXi ; tÞ  bðXi ; tÞdX ¼ 0 ðXi þ ui Þ  ½qðXi Þu

ðA:13Þ

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

32

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

References [1] Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397–434. [2] Camacho GT, Ortiz M. Computational modeling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899–938. [3] Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Numer Meth Eng 1999;44:1267–82. [4] Zhou F, Molinari JF. Dynamic crack propagation with cohesive elements: a methodology to address mesh dependence. Int J Numer Meth Eng 2004;59 (1):1–24. [5] Remmers JJC, De Borst R, Needleman A. A cohesive segments method for the simulation of crack growth. Comput Mech 2003;31:69–77. [6] Armero F, Garikipati K. An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids. Int J Solids Struct 1996;33(20–22):2863–85. [7] Belytschko T, Fish J, Englemann B. A finite element method with embedded localization zones. Comput Methods Appl Mech Eng 1988;70:59–89. [8] Samaniego E, Olive X, Huespe A. Contributions to the continuum modelling of strong discontinuities in two-dimensional solids. PhD thesis. International Center for Numerical Methods in Engineering, Monograph CIMNE; 2003 [No. 72]. [9] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45(5):601–20. [10] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46(1):133–50. [11] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int J Numer Meth Eng 2004;61(13):2316. 2316-2150. [12] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput Methods Appl Mech Eng 2007;196:2777–99. [13] Belytschko T, Lu YY. Element-free Galerkin methods for static and dynamic fracture. Int J Solids Struct 1995;32:2547–70. [14] Belytschko TT, Lu YY, Gu L. Crack propagation by element free Galerkin methods. Eng Fract Mech 1995;51(2):295–315. [15] Lu YY, Belytschko T, Tabbara M. Element-free Galerkin method for wave-propagation and dynamic fracture. Comput Methods Appl Mech Eng 1995;126(1–2):131–53. [16] Krysl P, Belytschko T. The element free Galerkin method for dynamic propagation of arbitrary 3-D cracks. Int J Numer Meth Eng 1999;44(6):767–800. [17] Ventura G, Xu J, Belytschko T. A vector level set method and new discontinuity approximations for crack growth by EFG. Int J Numer Meth Eng 2002;54 (6):923–44. [18] Bittencourt TN, Wawrzynek PA, Ingraffea AR, Sousa JL. Quasi-automatic simulation of crack propagation for 2D LEFM problems. Eng Fract Mech 1996;55(2):321–34. [19] Areias P, Belytschko T. Analysis of three-dimensional crack initiation and propagation using the extended finite element method. Int J Numer Meth Eng 2005;63(5):760–88. [20] Breitenfeld MS, Geubelle PH, Weckner O, Silling SA. Non-ordinary state-based peridynamic analysis of stationary crack problems. Comput Methods Appl Mech Eng 2014;272:233–50. [21] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 2000;48:175–209. [22] Foster JT, Silling SA, Chen WW. Viscoplasticity using peridynamics. Int J Numer Meth Eng 2010;81(10):1242–58. [23] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Comput Struct 2005;83:1526–35. [24] Silling SA, Epton M, Weckne O, Xu J, Askari E. Peridynamic states and constitutive modeling. J Elast 2007;88:151–84. [25] Silling SA, Lehoucq RB. Peridynamic theory of solid mechanics. Adv Appl Mech 2010;44:73–168. [26] Tupek MR, Rimodi JJ, Radovitzky R. An approach for incorporating classical continuum damage models in state-based peridynamics. Comput Methods Appl Mech Eng 2013;263:20–6. [27] Wang YT, Zhou XP, Xu X. Numerical simulation of propagation and coalescence of flaws in rock materials under compressive loads using the extended non-ordinary state-based peridynamics. Eng Fract Mech 2016;163:248–73. [28] Zhou XP, Wang YT, Qian QH. Numerical simulation of crack curving and branching in brittle materials under dynamic loads using the extended nonordinary state-based peridynamics. Eur J Mech A-Solid 2016;60:277–99. [29] Bobaru F, Yang MJ, Alves LF, Silling SA, Askari E, Xu JF. Convergence, adaptive refinement, and scaling in 1D peridynamics. Int J Numer Meth Eng 2009;77:852–77. [30] Silling SA. Dynamic fracture modeling with a meshfree peridynamic code. In: The second MIT conference on computational fluid and solid mechanics. MA: MIT; 2003. [31] Weckner O, Abeyaratne R. The effect of long-range forces on the dynamics of a bar. J Mech Phys Solids 2005;53:705–28. [32] Ha YD, Bobaru F. Studies of dynamic crack propagation and crack branching with peridynamics. Int J Fract 2010;162:229–44. [33] Agwai A, Guven I, Madenci E. Predicting crack propagation with peridynamics: a comparative study. Int J Fract 2011;171:65–78. [34] Ha YD, Bobaru F. Characteristics of dynamic brittle fracture captured with peridynamics. Eng Fract Mech 2011;78:1156–68. [35] Huang D, Lu GD, Qiao PZ. An improved peridynamic approach for quasi-static elastic deformation and brittle fracture analysis. Int J Mech Sci 2015;94– 95:111–22. [36] Huang D, Lu GD, Wang CW, Qiao ZP. An extended peridynamic approach for deformation and fracture analysis. Eng Fract Mech 2015;141:196–211. [37] Liu WY, Hong JH. Discretized peridynamics for linear elastic solids. Comput Mech 2012;50:579–90. [38] Prakash N, Seidel GD. A novel two-parameter linear elastic constitutive model for bond based peridynamics. In: 56th AIAA/ASCE/AHS/ASC structures, structural dynamics, and material conference, Kissimmee, Florida, USA; 2015. [39] Liu WY, Hong JH. Discretized peridynamics for brittle and ductile solids. Int J Numer Meth Eng 2012;89:1028–46. [40] Ha YD, Lee J, Hong JW. Fracturing patterns of rock-like materials in compression captured with peridynamics. Eng Fract Mech 2015;144:176–93. [41] Lee J, Liu WY, Hong JW. Impact fracture analysis enhanced by contact of peridynamic and finite element formulations. Int J Impact Eng 2016;87:108–19. [42] Lee J, Oh SE, Hong JW. Parallel programming of a peridynamics code coupled with finite element method. Int J Fract 2016. doi: http://dx.doi.org/ 10.1007/s10704-016-0121-y. [43] Stillinger FH, Weber TA. Computer simulation of local order in condensed phases of silicon. Phys Rev B 1985;31:5262–71. [44] Zhang ZN, Chen YX. Modeling nonlinear elastic solid with correlated lattice bond cell for dynamic fracture simulation. Comput Methods Appl Mech Eng 2014;279:325–47. [45] Pizzagalli L, Godet J, Guenole J, Brochard S, Holmstrom E, Nordlund K, et al. A new parametrization of the Stillinger-Weber potential for an improved description of defects and plasticity of silicon. J Phys Condens Matter 2013;25:055801. [46] Zhang ZN, Chen YX, Zheng H. A modified Stillinger-Weber potential-based hyperelastic constitutive model for nonlinear elasticity. Int J Solids Struct 2014;51:1542–54. [47] Zhang ZN, Ge XR. A new quasi-continuum constitutive model for crack growth in an isotropic solid. Eur J Mech A-Solid 2005;24:243–52. [48] Sau A. Peridynamic modeling of quasibrittle structure. The University of New Mexico; 2008. [49] Gerstle W, Sau N, Silling SA. Peridynamic modeling of concrete structures. Nucl Eng Des 2007;237:1250–8. [50] Griffiths DV, Mustoe GGW. Modelling of elastic continua using a grillage of structural elements based on discrete element concepts. Int J Numer Meth Eng 2001;50:1759–75. [51] Foster JT, Silling SA, Chen WN. An energy based failure criterion for use with peridynamic states. Int J Multiscale Comput 2011;9(6):675–87. [52] John R, Shah SP. Mixed-mode fracture of concrete subjected to impact loading. J Struct Eng – ASCE 1990;116:585–602. [53] John R. Mixed mode fracture of concrete subjected to impact loading [PhD dissertation]. Northwestern University; 1988.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031

X. Zhou et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

33

[54] Belytschko T, Tabbara M. Dynamic fracture using element-free Galerkin methods. Int J Numer Meth Eng 1996;39:923–38. [55] Song JH, Areias PMA, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. Int J Numer Meth Eng 2006;67:868–93. [56] Belytschko T, Chen H, Xu J, Zi G. Dynamic crack propagation based on loss of hyperbolicity with a new discontinuous enrichment. Int J Numer Meth Eng 2003;58:1873–905. [57] Gupta V, Rajagopal S, Gupta N. A comparative study of meshfree methods for fracture. Int J Damage Mech 2011;20:729–51. [58] Braun M, Fernández-Sáez J. A new 2D discrete model applied to dynamic crack propagation in brittle materials. Int J Solids Struct 2014;51:3787–97. [59] Ramulu M, Kobayashi AS. Mechanics of crack curving and branching – a dynamic fracture analysis. Int J Fract 1985;27:187–201. [60] Ravi-Chandar K. Dynamic fracture of nominally brittle materials. Int J Fract 1988;90:187–201. [61] Sharon E, Gross SP, Fineberg J. Local crack branching as a mechanism for instability in dynamic fracture. Phys Rev Lett 1995;74:5096–9. [62] Sharon E, Fineberg J. Microbranching instability and the dynamic fracture of brittle materials. Phys Rev B 1996;54:7128–39. [63] Fineberg J, Sharon E, Cohen G. Crack front waves in dynamic fracture. Int J Fract 2003;119:247–61. [64] Song J, Wang H, Belytschko T. A comparative study on finite element methods for dynamic fracture. Comput Mech 2008;45:239–50. [65] Wildman RA, Gazonas GA. A finite difference-augmented peridynamics method for reducing wave dispersion. Int J Fract 2014;90:39–52. [66] Zhou SJ, Lomdahl PS, Thomson R, Holian BL. Dynamic crack processes via molecular dynamics. Phys Rev Lett 1996;76:2318–21. [67] Abraham FF, Brodbeck D, Rafey RA, Rudge WE. Instability dynamics of fracture: a computer simulation investigation. Phys Rev Lett 1994;73:272–5. [68] Abraham FF, Brodbeck D, Rudge WE, Xu XP. A molecular dynamics investigation of rapid fracture mechanics. J Mech Phys Solids 1997;45:1595–619. [69] Fineberg J, Gross SP, Marder M, Swinney HL. Instability in dynamic fracture. Phys Rev Lett 1991;67:457–60. [70] Cramer T, Wanner A, Gumbsch P. Energy dissipation and path instabilities in dynamic fracture of silicon single crystals. Phys Rev Lett 2000;85:788–91. [71] Yang YF, Tang CA, Xia KW. Study on crack curving and branching mechanism in quasi-brittle materials under dynamic biaxial loading. Int J Fract 2012;177:53–72.

Please cite this article in press as: Zhou X et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.07.031