A micro-potential based Peridynamic method for deformation and fracturing in solids: A two-dimensional formulation

Available online at www.sciencedirect.com

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

A micro-potential based Peridynamic method for deformation and fracturing in solids: A two-dimensional formulation Jincheng Fana ,∗, Renwei Liub , Shaofan Lib , Xiurun Gea a

School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China b Department of Civil and Environmental Engineering, The University of California, Berkeley, CA 94720, USA Received 29 August 2019; received in revised form 8 November 2019; accepted 11 November 2019 Available online xxxx

Abstract A two-dimensional micro-potential based Peridynamic (MPPD) formulation is proposed to investigate deformation and fracturing in solid materials. To this end, square measure of the bond length change is applied and a generalized Peridynamic (PD) strain for plane stress and strain problems is developed to decompose the bond length change into two parts, one resulting from the volumetric PD strain and the other from the deviatoric PD strain. The micro-potential generated in the deformed bond is postulated as a function of both length change components. The nonlocal elastic strain energy density (NESED) at a material point is computed by the integral of the bond potential over the horizon. Through the Frechet derivative of the NESED, a general constitutive relation depending on the micro-potential function is well formulated. The model parameters are calibrated through the equivalence of the NESED with the linear elastic strain energy density. Several specific micro-potential functions and their corresponding constitutive force densities are carefully discussed. The discussions reveal that the proposed MPPD model not only can retrieve the well-known bond based PD model, but also allows specific materials with any Poisson ratio to be represented. Moreover, one failure criterion based on the total micro-bond potential is developed for the present model. The reliability and efficiency of the MPPD model is demonstrated through numerical tests. Simulation results show that the proposed model is capable of investigating deformation and cracking in solids. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Peridynamics; Micro-potential; Constitutive model; Fracturing; Poisson ratio

1. Introduction Peridynamics (PD) is a new nonlocal formulation of solid mechanics in terms of a space integral and time differential equation instead of the space differential way used in the classical theory of solid mechanics [1–3]. The integral equation makes sense even in the discontinuities such as the cracks and voids, which makes PD a promising nonlocal theory that could unify discontinuities and continuities in a framework and a suitable method for studying the complex fracturing in solid materials without treating defects. It has been shown in [4,5] that the PD equation converges to the Navier equation of linear elasticity as the nonlocal horizon goes to zero in the limit. Besides, Seleson et al. [6] revealed that PD could play as an upscaling of molecular mechanics in some case. ∗ Corresponding author.

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

2

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

Practical applications of PD have been demonstrated in the study of a variety of materials such as composites [7–9], crystals [10,11], glass [12,13], geomaterials [14–17], to name a few. In the original bond based Peridynamics (BPD), the mechanical response of each bond depends only on the deformation of itself and one single micro-elastic parameter is used to characterize the linear force–stretch relationship. As a result, the applicability of the original BPD is limited to isotropic materials with the fixed Poisson ratio of 1/4 for 3D and plane strain problems and 1/3 for plane stress problems [18]. To effectively apply PD to study the failure behavior of more materials with Poisson ratio other than 1/3 or 1/4, several attempts have been made [2,19–22]. Wang et al. [20] introduced a novel conjugated bond-pair based PD model in which the bond force was not only related to the stretch but also related to the rotation of conjugated bonds. Hu and Madenci [21] classified the bonds as the normal bonds and shear bonds. The bond parameters were determined by the elastic constants through equating the PD force with that obtained by the classical continuum mechanics under the same deformation. The method proposed by Hu and Madenci [21] has been used to model the deformation and failure of composite laminates with broader Poisson ratios. In fact, Silling et al. [2] developed another version of PD referred to as the ordinary state based PD (OSPD) to solve the fixed Poisson ratio problem. A nonlocal elastic strain energy density (NESED) function expressed in terms of bulk modulus and shear modulus was postulated. The bond force density was obtained through the Frechet derivative of the NESED with respect to the deformation state. Ren et al. [22] developed a distinct formulation for the NESED and derived the constitutive relation by the derivative of the NESED with respect to the isotropic and deviatoric parts of the displacement. In the present study, a micro-potential based PD (MPPD) constitutive model is proposed for deformation and fracturing in solids. The paper is organized as follows: in the next section, the basic theory of PD is briefly reviewed with special regard to the BPD and OSPD. In Section 3, a general two-dimensional micro-potential based PD constitutive formulation is derived. Several specific potential functions as well as the corresponding constitutive relations are introduced and analytically discussed. In Section 4, a failure criterion based on the total micro-bond potential is proposed for the MPPD model. Numerical validation of the present model is performed in Section 5. At last, we present the concluding remarks. 2. A brief review of PD In general, PD can be classified as the BPD and the state based Peridynamics (SPD) according to whether the mechanical response of a material point in the solid depends on the collective deformation of all bonds in the horizon or not [2,3]. 2.1. BPD In the BPD, a material point x is assumed to interact with any other material point x′ that lies in a finite domain around x. The interaction between the pair points is called as the pair force density. The equation of motion at the material point x is formulated in the form of being integral in space and differential in time: ∫ ρ u¨ (x, t) = f (u(x′ , t) − u(x, t), x′ − x) d Vx′ + b(x, t), (1) H

where ρ is the mass density, u is the displacement vector field, u¨ denotes the second derivative of the displacement, b is the }prescribed body force and f represents the pair wise force density. H is defined as H = { ′ x ∈ B||x′ − x| < δ , in which B denotes the material body and δ is a positive constant referred to as the horizon. H represents a circular domain centered at x with a radius of δ for 2D and a spherical domain with a radius of δ for 3D. Note that the pair force density between material points with the distance beyond δ is assumed to be zero in PD. The conservation of linear and angular momentum places two admissible conditions on the pair force density: f (−η, −ξ ) = − f (η, ξ ) ,

(2)

(ξ + η) × f (η, ξ ) = 0,

(3)

where (ξ =) x′ − x denotes the relative position of two material points in the reference configuration, and η = u x′ , t − u (x, t) represents the relative displacement of the two points. In PD, ξ is called the initial bond.

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

3

For an isotropic elastic medium, the pairwise force density f is determined as: f (η, ξ ) =

∂w(η, ξ ) , ∂η

(4)

where w (η, ξ ) denotes the micro-potential, which characterizes the micro-elastic energy density per unit volume stored in the deformed bond. In the BPD, the micro-potential is defined as: 1 w = cs 2 ξ, (5) 2 in which c is the micro-elastic parameter and s is the bond stretch: |ξ + η| − |ξ | . |ξ | Inserting Eq. (5) into Eq. (4) leads to: s=

(6)

ξ +η . (7) |ξ + η| Eq. (7) shows that the pair force density in the original BPD depends linearly on the bond stretch and the direction of the force density is collinear with that of the deformed bond. Obviously, such pair force density satisfies the admissibility conditions. The NESED at a material point is computed by the integral of the micro-potential stored in the bonds connected to that point: ∫ 1 W = w(η, ξ )d Vξ . (8) 2 H Inserting Eq. (5) into Eq. (8) and equating the obtained NESED with the linear elastic strain energy density (LESED), the micro-elastic constant c is calibrated as [18]: ⎧ 18k 3D ⎪ ⎪ ⎨ πδ4 9k plane stress (9) c= πδ 3 ⎪ ⎪ ⎩ 72k plane strain. f (η, ξ ) = cs

5πδ 3

Note that there are two independent elastic constants for isotropic materials in the classical continuum mechanics. However, the micro-potential in Eq. (5) contains only one microscale parameter c, and thus the computed NESED contains one single parameter as well. Thus, the BPD is limited to materials with the fixed Poisson ratio of 1/4 for 3D and plane strain problems and 1/3 for plane stress problems. In [23], damage is introduced through the elimination of the interaction between two material points if the critical bond stretch is reached. Crack surface is formed when all interactions across one plane is eliminated. By equating the total work required to eliminate all interactions across the crack surface with the critical energy release rate, the critical stretch is analytically expressed as [18]: ⎧ √ 5G 0 ⎪ ⎪ 3D ⎪ 9kδ ⎪ ⎨ √ 4π G 0 s0 = (10) plane stress 9kδ ⎪ ⎪ √ ⎪ ⎪ ⎩ 5π G 0 plane strain. 18kδ Local damage at a material point x is defined as: ∫ Fail(x, x′ , t)d Vξ φ(x, t) = 1 − H ∫ , H d Vξ

(11)

where Fail(x, x′ , t) =

{

1 0

|x − x′ | < s0 for all t ′ ≤ t else.

(12)

4

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

2.2. OSPD The notion of “state” is introduced in the SPD to describe the collectively determined quantities. A state A⟨·⟩ is defined as a map function on H. The angle bracket encloses the bond vector. A is called a vector state if the value of A⟨·⟩ is a vector. If the operation value is a scalar, it is referred to as a scalar state. In the SPD, the deformation state Y is introduced to represent the collective deformation of the bonds connected to the same material point. Under its operation, any bond ξ in the reference configuration is mapped to the deformed bond Y ⟨ξ ⟩ in the current configuration. Accordingly, another notation called the force state T is defined to characterize the forces density within the bonds. The relation between the force state and the deformation state constitutes the constitutive model in the SPD framework. According to whether the direction of the bond force density is parallel to the deformed bond or not, the SPD can be classified as the ordinary SPD (OSPD) and the non-ordinary SPD (NOSPD). Here we only give a brief review of the OSPD and one can refer to [2,3] for the introduction of the NOSPD. The fundamental equation of motion is expressed as: ∫ { ⟨ ⟩ ⟨ ⟩} ρ(x)u¨ (x) = T [x] x′ − x − T [x′ ] x − x′ d Vx′ + b(x). (13) H

In the framework of the OSPD, the nonlocal strain energy density W is a functional of deformation state and written as: W = Wˆ (Y ) .

(14)

The constitutive relation is defined as the Frechet derivative of Wˆ : T = Wˆ Y ,

(15)

where the Frechet derivative is introduced as: Wˆ (Y + δY ) − Wˆ (Y ) = Wˆ Y • δY + o (∥δY ∥) ,

(16) √

in which δY denotes an increment in the deformation state and ∥δY ∥ = δY • δY . In the original OSPD model for plane stress problems [24], a notion of nonlocal dilatation is introduced as: θ=

2(1 − 2ν) (ωx) • e , 1 − ν (ωx) • x

(17)

in which ω is a scalar influence state and e is the extension state: e ⟨ξ ⟩ = |Y ⟨ξ ⟩| − x ⟨ξ ⟩ = |Y ⟨ξ ⟩| − |ξ | .

(18)

The NESED for plane stress problems is defined as: W (Y ) =

k′ 2 α θ + (ωed ) • ed , 2 2

µ(1+ν) in which k ′ = k + 9(1−2ν) 2, α = deviatoric part of extension state: 2

(19) 8µ (ωx)•x

where k and µ are the bulk and shear modulus respectively, ed is the

θ |ξ |. 3 The derivative of Eq. (19) yields the following constitutive relation: [ ] ) ω⟨ξ ⟩|ξ | 2(1 − 2ν) ( ′ α Y ⟨ξ ⟩ d d T ⟨ξ ⟩ = + αω⟨ξ ⟩e ⟨ξ ⟩ . k θ − (ωe ) • x 1−ν 3 (ωx) • x |Y ⟨ξ ⟩| ed ⟨ξ ⟩ = e ⟨ξ ⟩ − ei ⟨ξ ⟩ ,

ei ⟨ξ ⟩ =

(20)

(21)

It is readily to prove that the constitutive relation in Eq. (21) satisfies the conservations of linear and angular momentums. One can refer to [2] for the 3D OSPD model and [25] for an axisymmetric OSPD model. In the original BPD and OSPD, the horizon size is set constant. For varying horizon size, one can refer to the dual-horizon model introduced in [26,27].

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

5

3. Proposed constitutive model In the original BPD, micro-potential with the magnitude of 1/2cs 2 ξ is generated if one bond has a stretch of s, see Eq. (5). The micro-potential in the bond connected to the material point x contributes to the macroscopic NESED, see Eq. (8). The pair force density is computed by the derivative of the micro-bond potential with respect to the relative displacement, see Eq. (4). The derived pair force density depends linearly on the bond stretch like an elastic “spring”, see Eq. (7). It can be considered that the micro-potential could reveal and determine the microscopic deformation mechanism. However, the simple micro-potential function with a single material constant leads to the fixed Poisson ratio problem to some extent. In [24], the modeling strategy is: firstly, one volumetric strain formulation, e.g. Eq. (17) is proposed. Secondly, a NESED function, e.g. Eq. (19) is postulated from the macroscopic point of view. Thirdly, the constitutive relation is obtained by the Frechet derivative of the NESED. Since the postulated NESED contains two material constants, the fixed Poisson ratio problem confronted by the original BPD could be fixed. However, the notion of micro-bond potential is abandoned in the OSPD introduced in [24]. In other words, the origin of the macroscopic NESED is no longer explicitly explained. Combine the modeling advantages of the BPD and OSPD, we propose a novel micro-potential based PD formulation in the OSPD framework for deformation and fracturing in solids. 3.1. General constitutive relation As the solid deforms, the bond length may change. In the present analysis, square measure of the bond length change is applied: ∆ = |Y ⟨ξ ⟩|2 − |ξ |2 ,

(22)

where Y ⟨ξ ⟩ and ξ denotes the bond in the current and reference configuration respectively. In plane stress and strain problems, the initial bond is expressed as ξ = {ξ1 , ξ2 , 0} with ξ1 = ξ cos(θ ) and ξ2 = ξ sin(θ ). According to the Cauchy–Born rule, the deformed bond can be expressed as: Y ⟨ξ ⟩ = Fξ ,

(23)

where F is the deformation gradient. Substituting Eq. (23) into Eq. (22) yields: ∆ = 2ξ T Eξ ,

(24) 1 2

(

T

)

F F−1 . in which E is the Green strain tensor and E = For plane stress and strain problems, the Green strain can be written in the matrix form as: ⎤ ⎤ ⎤ ⎡ ⎡ ⎡ e11 e12 0 ε11 ε12 0 1 0 0 θ [E] = ⎣ ε21 ε22 0 ⎦ = ⎣ 0 1 0 ⎦ + ⎣ e21 e22 0 ⎦ , 3 0 0 1 0 0 e33 0 0 ε33

(25)

ν in which ei j = εi j − 13 θ δi j (i, j = {1, 2}), for plane stress problems: ε33 = − 1−ν (ε11 + ε22 ), θ = 1−2ν (ε11 + ε22 ), 1−ν 1 1+ν e33 = − 3(1−ν) (ε11 + ε22 ), and for plane strain problems: ε33 = 0, θ = ε11 + ε22 , e33 = − 3 (ε11 + ε22 ). Inserting Eq. (25) into Eq. (24), we can get:

∆ = ∆v + ∆d ,

(26)

in which 2θ 2 ξ , ∆d = 2ei j ξi ξ j . (27) 3 ∆v and ∆d are considered as the length change components resulting from the volume deformation and deviatoric deformation respectively. In the framework of the SPD, the sole quantity to describe the deformation in solids is the deformation state Y . In order to compute the bond length change components, we propose the formulation of generalized PD strain components as: ( )( ) ∫ ξi ξ j 1 Y ⟨ξ ⟩ · Y ⟨ξ ⟩ 1 εi j = ω⟨ξ ⟩ −1 2 − δi j d Sξ , (28) A H ξ ·ξ ξ ·ξ 2 ∆v =

6

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

∫ in which i, j = {1, 2}, A = H ω⟨ξ ⟩d Sξ and d Sξ = ξ dξ dθ . It is noted that A is introduced to delete the surface effect. The equivalence of the proposed PD strain components with the classical ones is proven in Appendix A. According to Eq. (28), the volumetric PD strain is computed as: ( ) ∫ B Y ⟨ξ ⟩ · Y ⟨ξ ⟩ θ= ω⟨ξ ⟩ − 1 d Sξ , (29) A H ξ ·ξ for plane stress problems and B = 1 for plane strain problems. in which B = 1−2ν 1−ν The deviatoric PD strain components are formulated as: ( )( ) ∫ ξi ξ j 1 Y ⟨ξ ⟩ · Y ⟨ξ ⟩ 3 + 2B ei j = ω⟨ξ ⟩ −1 2 − δi j d Sξ . A H ξ ·ξ ξ ·ξ 6

(30)

In the original BPD, when the bond length changes, micro-potential is generated and contributes to the macroscopic NESED at that material point, see Eq. (8). In the present analysis, the square measured bond length change is decomposed into two parts, one resulting from the volumetric PD strain and the other from the deviatoric PD strain. Remember that the LESED in the classical continuum mechanics is expressed as the function of volumetric strain and deviatoric strain. As a result, we postulate that the micro-potential in Y ⟨ξ ⟩ is owned by the material point x and takes the following form: w(x, ξ , Y ) = w(∆ ˆ v , ∆d ),

(31)

in which ∆v = ∆v (θ ) and ∆d = ∆d (ei j ). Remark. In the original BPD, the micro-potential in the deformed bond is just related to the initial bond length and stretch of itself. The stored micro-potential is shared by the two material points connected by that bond. In the present analysis, the micro-potential owned by x not only depends on the initial bond length, but also is related to its direction in the reference configuration. Moreover, the proposed micro-potential is a functional of the deformation state Y , which indicates that the individual bond potential tends to be affected by the collective deformation in the horizon. Specific micro-potential functions will be discussed later. The NESED at x is computed by the integral of w over the horizon: ∫ w(x, ξ , Y )d Vξ . W =

(32)

H

It should be noted that the coefficient of 1/2 in Eq. (8) is omitted here, since w defined in the present analysis represents the micro-bond potential in Y ⟨ξ ⟩ that is solely owned by the material point x. The variation of the NESED in Eq. (32) is expressed as: ) ∫ ∫ ( ∂w ∂w δwd Vξ = t δW = δ∆v + δ∆d d Sξ , (33) ∂∆d H H ∂∆v in which t represents the thickness and is unit of one. According to Eq. (26), we can obtain: δ∆d = δ∆ − δ∆v . Substituting the above equation into Eq. (33) yields: ) ] ∫ [( ∂w ∂w ∂w δW = t − δ∆v + δ∆ d Sξ . ∂∆v ∂∆d ∂∆d H According to Eqs. (27) and (29), the variation of ∆v can be calculated as: ∫ 2ξ 2 4B 2 Y ⟨β⟩ · δY ⟨β⟩ δ∆v = δθ = ξ ω⟨β⟩ d Sβ . 3 3A β ·β H

(34)

(35)

(36)

Due to Eq. (22), the variation of the total length change is expressed as: δ∆ = 2Y ⟨ξ ⟩ · δY ⟨ξ ⟩.

(37)

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

Inserting Eqs. (36) and (37) into Eq. (35), we can have: ) ∫ ( ∫ ∫ 4B ∂w Y ⟨β⟩ · δY ⟨β⟩ ∂w ∂w δW = − ξ 2 d Vξ ω⟨β⟩ d Sβ + 2 Y ⟨ξ ⟩ · δY ⟨ξ ⟩d Vξ . 3A H ∂∆v ∂∆d β · β ∂∆ d H H

7

(38)

Exchanging the dummy integration variables, ξ ↔ β in the first term on the right hand, Eq. (38) can be rewritten as: δW =

∫

[

H

4B ω⟨ξ ⟩ξ −2 3A

∫

(

H

∂w ∂w − ∂∆v ∂∆d

)

β 2 d Sβ + 2

] ∂w Y ⟨ξ ⟩ · δY ⟨ξ ⟩d Vξ . ∂∆d

According to Eq. (15), the proposed constitutive relation is formulated as: [ ) ] ∫ ( 4B ∂w ∂w ∂w T ⟨ξ ⟩ = ω⟨ξ ⟩ξ −2 − β 2 d Sβ + 2 (ξ ) Y ⟨ξ ⟩. 3A ∂∆d ∂∆d H ∂∆v

(39)

(40)

The bond force density in Eq. (40) can be decomposed as: T ⟨ξ ⟩ = T 1 ⟨ξ ⟩ + T 2 ⟨ξ ⟩,

(41)

where [∫ ( ) ] ∂w ∂w 4B −2 2 ω⟨ξ ⟩ξ − β d Sβ Y ⟨ξ ⟩, T 1 ⟨ξ ⟩ = 3A ∂∆d H ∂∆v ∂w (ξ )Y ⟨ξ ⟩. T 2 ⟨ξ ⟩ = 2 ∂∆d

(42) (43)

Remark. The integral in Eq. (42) represents the contribution of the collective responses in the horizon to the individual bond force density. Besides, the direction of T 1 and T 2 are parallel to that of Y ⟨ξ ⟩, which indicates that the conservations of linear and angular momentums are satisfied. In general, the integral in Eq. (42) will not be ∂w ∂w ≡ ∂∆ identically equal to zero unless ∂∆ , and thus the proposed MPPD formulation can be considered as a kind v d ∂w ∂w of OSPD models. As ∂∆v ≡ ∂∆d , the proposed MPPD model can be reduced to a kind of BPD models for some specific micro-potential function such as that given in Eq. (44). 3.2. Specific constitutive relation-1 We consider the micro-potential in Eq. (31) as a quadratic function of ∆: wˆ = ω1 ⟨ξ ⟩(∆v + ∆d )2 = ω1 ⟨ξ ⟩∆2 .

(44)

The partial derivatives of wˆ with respect to ∆v and ∆d are written as: ∂ wˆ ∂ wˆ = = 2ω1 ⟨ξ ⟩∆. (45) ∂∆v ∂∆d Inserting Eq. (45) into Eqs. (42) and (43) gives the bond force density corresponding to the micro-potential in Eq. (44) as: T 1 ⟨ξ ⟩ = 0,

T 2 ⟨ξ ⟩ = 4∆ω1 ⟨ξ ⟩Y ⟨ξ ⟩.

(46)

The equivalence of the NESED with the LESED under homogeneous deformation is applied to calibrate ω1 ⟨ξ ⟩. In the case of homogeneous deformation, Eq. (44) can be expressed as: wˆ = 4ω1 ⟨ξ ⟩εi j ξi ξ j εkl ξk ξl , and the NESED is formulated as: (∫ ) W PD = ω1 ⟨ξ ⟩ξ 4 d Vξ L i jkl εi j εkl ) (∫H (∫ ) ( 2 ) 1 4 2 2 4 ω1 ⟨ξ ⟩ξ d Vξ (ε11 + ε22 )2 , = ω1 ⟨ξ ⟩ξ d Vξ ε11 + ε22 + 2ε12 + 2 H H in which L i jkl is given in Appendix B.

(47)

(48)

8

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

The LESED function for plane stress and strain problems is formulated as: ) 3k(1 − 2ν) ν ′ 3k(1 − 2ν) ( 2 2 2 + W classical = ε11 + ε22 + 2ε12 (ε11 + ε22 )2 , 2(1 + ν) 2(1 + ν) 1 − ν ′ in which ν ′ = ν for plane stress problems and ν ′ = ν/(1 − ν) for plane strain problems. According to Eqs. (48) and (49), we can obtain: { 3k { 1 ∫ plane stress plane stress 8 3 ω1 ⟨ξ ⟩ξ 4 d Vξ = , ν= . 3k 1 plane strain plane strain H 5 4

(49)

(50)

Remark. The model corresponding to the micro-potential in Eq. (44) is only applicable for materials with the fixed Poisson ratio, ν ≡ 1/3 for plane stress problems and ν ≡ 1/4 for plane strain problems. For plane stress problems, to satisfy Eq. (50) ω1 is written as: ω1 ⟨ξ ⟩ =

3k ω3 ⟨ξ ⟩ξ −2 . 8(ω3 x) • x

(51)

In the following, four distinct functions ω3 will be introduced and their corresponding pair force density functions are formulated. • Case A: ω3 ⟨ξ ⟩ = ξ −1 .

(52)

Substituting Eqs. (52) and (51) into Eq. (46) leads to: f (ξ ) = 2T 2 ⟨ξ ⟩ =

9k Y ⟨ξ ⟩ s(1 + s)(2 + s) . 2π tδ 3 |Y ⟨ξ ⟩|

(53)

As the deformation is small, |Y ⟨ξ ⟩| ≈ ξ . The above bond force density can be approximated as: f (ξ ) =

9k Y ⟨ξ ⟩ s , π δ 3 |Y ⟨ξ ⟩|

(54)

which is consistent with the original BPD pair force formulation. Remark. In this case, the obtained bond force density is just related to the bond stretch and is not concerned with the initial bond length. Besides, the BPD model developed by Silling [1] can be approximately retrieved by the present MPPD model under small deformation. • Case B: ω3 ⟨ξ ⟩ ≡ constant, which results in: ω3 ⟨ξ ⟩ 2 = . (ω3 x) • x π tδ 4

(55)

(56)

Substituting Eqs. (56) and (51) into Eq. (46) results in: f (ξ ) = 2T 2 ⟨ξ ⟩ =

6k ξ Y ⟨ξ ⟩ . s(1 + s)(2 + s) π tδ 3 δ |Y ⟨ξ ⟩|

(57)

• Case C: ξ2 ). δ2 Inserting Eqs. (58) and (51) into Eq. (46) yields: ω3 ⟨ξ ⟩ = exp(−

f (ξ ) = 2T 2 ⟨ξ ⟩ =

3k ξ ξ 2 Y ⟨ξ ⟩ s(1 + s)(2 + s) exp(− 2 ) . −1 3 (1 − 2e )π tδ δ δ |Y ⟨ξ ⟩|

(58)

(59)

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

9

Fig. 1. The variation of the pair force density for different initial bond lengths.

• Case D: 10 (δ + 3ξ )(δ − ξ )3 . π δ6 Substituting Eqs. (60) and (51) into Eq. (46) gives: ω3 ⟨ξ ⟩ =

(60)

84k ξ ξ ξ Y ⟨ξ ⟩ s(1 + s)(2 + s) (1 + 3 )(1 − )3 . (61) 3 π tδ δ δ δ |Y ⟨ξ ⟩| The variations of four pair force densities for different initial bond lengths are illustrated in Fig. 1. The pair force density in case A does not depend on the initial bond length. The pair force density in case B exhibits a linear dependence on the initial bond length. Both pair force densities in case C and case D increase first and then decrease as the initial bond length increases. Under the same stretch, the pair force density of the bond with the largest length reaches the maximum in case B, is slightly smaller than the peak in case C and approaches zero in case D. f (ξ ) = 2T 2 ⟨ξ ⟩ =

3.3. Specific constitutive relation-2 Introducing the following quadratic micro-potential function: wˆ = ω1 ⟨ξ ⟩∆2v + ω2 ⟨ξ ⟩∆2d . To satisfy the equivalence of the NESED with the LESED, ω1 ⟨ξ ⟩ and ω2 ⟨ξ ⟩ should obey: { µ(1+ν)(7−11ν) ∫ ∫ plane stress 8(1−2ν)2 4 ω1 ⟨ξ ⟩ξ d Vξ = =C , ω2 ⟨ξ ⟩ξ 4 d Vξ = µ. µ(7+4ν) H H plane strain 8(1−2ν)

(62)

(63)

One can refer to Appendix B for the proof. The partial derivatives of wˆ with respect to ∆v and ∆d are computed as: ∂ wˆ 4θ 2 = 2ω1 ⟨ξ ⟩∆v = ξ ω1 ⟨ξ ⟩, ∂∆v 3 ∂ wˆ = 2ω2 ⟨ξ ⟩∆d = 4ω2 ⟨ξ ⟩ei j ξi ξ j . ∂∆d To satisfy the necessary condition given in Eq. (63), we have: C µ ω1 ⟨ξ ⟩ = ω3 ⟨ξ ⟩ξ −2 , ω2 ⟨ξ ⟩ = ω4 ⟨ξ ⟩ξ −2 . (ω3 x) • x (ω4 x) • x

(64) (65)

(66)

10

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

The constitutive relation corresponding to the micro-potential expressed in Eq. (62) can be obtained by substituting the above equations into Eqs. (42)–(43). 3.4. Specific constitutive relation-3 For geomaterials such as rocks, concrete and ceramics, the failure process is dominated by cohesive fracturing. To model the cohesive deformation, we consider the following exponential–exponential micro-potential function: ( ) wˆ = φ − φ (1 + ω1 ⟨ξ ⟩|∆v |) exp (−ω1 ⟨ξ ⟩|∆v |) exp −(ω2 ⟨ξ ⟩∆d )2 , (67) in which φ is a material constant with the dimension of MT−2 L−4 , ω1 ⟨ξ ⟩ and ω2 ⟨ξ ⟩ are two functions that need to satisfy the following conditions: { µ(1+ν)(7−11ν) ∫ ∫ plane stress 8(1−2ν)2 2 4 , φ ω22 ⟨ξ ⟩ξ 4 d Vξ = µ, (68) φ ω1 ⟨ξ ⟩ξ d Vξ = µ(7+4ν) H H plane strain 8(1−2ν) by which the equivalence of the NESED with the LESED under infinitesimal deformation can be reached, see Appendix B. The partial derivatives of the micro-potential function (67) with respect to ∆v and ∆d are derived as: ( ) ∂ wˆ = φω12 ⟨ξ ⟩∆v exp (−ω1 ⟨ξ ⟩|∆v |) exp −(ω2 ⟨ξ ⟩∆d )2 , (69) ∂∆v ( ) ∂ wˆ = 2φω22 ⟨ξ ⟩∆d (1 + ω1 ⟨ξ ⟩|∆v |) exp (−ω1 ⟨ξ ⟩|∆v |) exp −(ω2 ⟨ξ ⟩∆d )2 . (70) ∂∆d The constitutive relation corresponding to the exponential–exponential micro-potential can be obtained by inserting Eqs. (69) and (70) into Eqs. (42) and (43). Due to the complex exponential function, it is difficult to get an explicit formulation for the bond force density in most cases. Here we discuss the simple case that the solid deforms homogeneously and isotropically. The plane stress problem is considered, i.e. ε11 (x, y) = ε22 (x, y) = ε. According to Eq. (27), the bond length changes are computed as: 2(1 + ν)ε 2 4(1 − 2ν)ε 2 ξ , ∆d = ξ . 3(1 − ν) 3(1 − ν) For illustration, we choose: 1 1 ∫ ∫ ω1 ⟨ξ ⟩ = , ω2 ⟨ξ ⟩ = , γ εc ξ 2 t H d Sξ ζ εc ξ 2 t H d Sξ ∆v =

(71)

(72)

in which εc is the limit uniaxial compressive strain, γ and ζ are two dimensionless constants that can be determined by Eq. (68) as: φ µ(1 + ν)(7 − 11ν)εc2 = , 2 γ 8(1 − 2ν)2

φ = µεc2 . ζ2

(73)

Following Eq. (73), we can obtain: (1 + ν)(7 − 11ν) ζ2 = . 2 γ 8(1 − 2ν)2

(74)

Inserting Eqs. (71)–(74) into Eqs. (69) and (70) leads to: ∂ wˆ µ(1 + ν)(7 − 11ν) 1 ∫ = Ψ1 (ε), ∂∆v 6(1 − ν)(1 − 2ν)t H d Sξ ξ 2 ∂ wˆ 4µ(1 + ν) 1 ∫ = Ψ2 (ε), ∂∆d 3(1 − ν)t H d Sξ ξ 2

(75) (76)

in which ( ) 4(1 − 2ν) 4(1 + ν)2 2 , Ψ1 (ε) = ε exp − |ε| − 2 ε 3γ (1 − ν)εc 9ζ (1 − ν)2 εc2

(77)

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

[ ] 4(1 − 2ν) Ψ2 (ε) = 1 + |ε| Ψ1 (ε). 3γ (1 − ν)εc Substituting Eqs. (75) and (76) into the integral term in Eq. (42) yields: ) ∫ ( ∂ wˆ ∂ wˆ µ(1 + ν) − β 2 d Sβ = Ψ (ε), ∂∆d 6(1 − ν)t H ∂∆v

11

(78)

(79)

in which [ Ψ (ε) =

] 32(1 − 2ν) 5ν − 1 − |ε| Ψ1 (ε). 1 − 2ν 3γ (1 − ν)εc

Thus, the bond force density resulting from the isotropic compression can be computed as: 2µ(1 − 2ν)(1 + ν) |T 1 ⟨ξ ⟩| = ω⟨ξ ⟩ξ −1 (1 + ε)Ψ (ε), 9A(1 − ν)2 t (∫ )−1 8µ(1 + ν) d Sξ ξ −1 (1 + ε)Ψ2 (ε). |T 2 ⟨ξ ⟩| = 3(1 − ν)t H

(80)

(81) (82)

The bond force density calculated by other available models are expressed as: • the bond force density computed by Le et al. [24] is: |T ⟨ξ ⟩| =

4µ(1 + ν)ω⟨ξ ⟩ξ (∫ ) ε. (1 − ν)t H ω⟨ξ ⟩ξ 2 d Sξ

(83)

• the bond force density formulated by Madenci [28] is: |T ⟨ξ ⟩| =

4µ(1 + ν − 4ν 2 )ω⟨ξ ⟩ξ (∫ ) ε. (1 − ν)2 t H ω⟨ξ ⟩ξ 2 d Sξ

(84)

• the bond force density determined by Gerstle et al. [18] is: 9E ε. (85) |T ⟨ξ ⟩| = 2tπ δ 3 It should be noticed that the original BPD model is restricted to materials with the fixed Poisson ratio of 1/3 for plane stress problems. Thus, Eq. (85) is only applicable for materials with Poisson ratio of 1/3. For comparison, the bond force density computed by our model and the other models, i.e. Eqs. (81)–(85) are discussed in the case 2 that ν = 1/3. The influence functions in all models are chosen as: ω⟨ξ ⟩ = exp(− ξδ2 ). As a result, the bond force density computed by the proposed model, i.e. Eqs. (81) and (82) is: [ ( 2) ] |T ⟨ξ ⟩| 1 ξ δ = exp − (1 + ε). (86) Ψ (ε) + 2Ψ (ε) 2 3 2 E/(tπ δ ) 12(1 − exp(−1)) δ ξ The bond force density determined by Eqs. (83) and (84) reads: ( 2) |T ⟨ξ ⟩| 3 ξ ξ = exp − 2 ε. E/(tπ δ 3 ) 1 − 2 exp(−1) δ δ

(87)

The bond force densities computed by Eqs. (85) and (87) increase linearly with the increasing |ε|. However, a distinct relation is exhibited by the proposed model. For illustration, we take the bond with the initial length of ξ = δ/3 for example. εc = 0.003 is assumed. The bond force density determined by the proposed model is depicted in Fig. 2. At the first stage, the force density increases linearly with a slope close to that computed by Gerstle et al.’s model. Then, the bond force density goes nonlinearly to the peak, followed by a gradual reduction to zero. The maximum is close to the force density by Le et al.’s model at the stationary strain of the proposed model. As ε = 0.002, the variations of the bond force densities in the horizon are illustrated in Fig. 3. The bond force density computed by Gerstle et al.’s model remains constant in the horizon. The bond force density calculated by Le et al.’s model is related to the initial bond length. Although all bonds in the horizon deform with the identical strain, the resulted bond force densities are not identical. As the initial bond length increases, the force density increases slowly and reaches the peak at ξ = 0.7δ, after which it begins to decrease. For the proposed model, the bond force density tends to be infinite as the initial bond length goes to zero. The force density of the bond with the initial

12

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

Fig. 2. The variation of the bond force density for different |ε|.

Fig. 3. The variation of the bond force density for different initial bond lengths.

length of ξ < 0.2δ decreases sharply as the length increases. If the bond length further increases, the force density decreases gradually and finally approaches a limit. It concludes that under identical deformation, the material point away from the studied point applies smaller force density. 4. Failure criterion In this section, the failure criterion for the MPPD model will be introduced. As shown in Fig. 4, the initial body Ω0 is mapped to the deformed body Ωn by the deformation state Y n . Due to deformaiton, micro-potential is generated and stored in the deformed bond that connects x and x′ . As pointed above, the micro-potential owned by the material point x is computed by w(x, ξ , Y n ) and the micro-potential owned by x′ is w(x′ , −ξ , Y n ). As a result, the total micro-potential stored in the deformed bond is computed as: w(x, x′ , Y n ) = w(x, ξ , Y n ) + w(x′ , −ξ , Y n ).

(88)

The proposed failure criterion states that as the total micro-potential exceeds the critical value wc , the bond connecting x and x′ breaks. The failing function is expressed as: { 1 w(x, x′ , Y n ) < wc for all t ′ ≤ t ′ Fail(x, x , t) = (89) 0 else.

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

13

Fig. 4. Deformation map.

5. Numerical applications In Section 3, the general MPPD formulation has been developed, and several specific micro-potential functions as well as the corresponding constitutive relations have been analytically discussed. In this section, numerical tests will be conducted to validate the efficiency of the MPPD for deformation and fracturing in solid materials. 5.1. Square plate with a hole The first test is to compute the deformation in an elastic square plate with a centered circular hole, see Fig. 5. The length L of the plate is 1 m. The density and elastic properties of the plate are taken as: ρ = 2265 kg/m3 , E = 30 GPa and ν = 0.2. Tensile stress with the magnitude of 200 MPa is uniformly applied on vertical boundaries. Different hole sizes are studied. The hole diameter of Plate-1 is 0.2 m, i.e. L/D = 5. The hole diameter of Plate-2 is 0.05 m, i.e. L/D = 20. The plane stress model proposed in Section 3.3 is applied to investigate the deformation in the plate. Methods have been developed to obtain the solutions for PD problems subject to static or quasi-static conditions, see [29,30]. The present analysis is conducted using the adaptive dynamic method introduced in [29] with a total of 5000 steps. The time step is ∆t = 1 s. Plate-1 is discretized by uniformly distributed 501 × 501 particles and the space between two neighboring particles is 0.002 m. The horizon size δ is set to be 0.006 m. The velocities and displacements at two specific points, i.e. point A (L/4, 0) and point B (0, L/4) are recorded with computational steps, see Fig. 6, from which we can get that the computation is well converged after 2000 steps. Velocities approach zero at the step of 2000 and afterwards. The horizontal displacement at point A converges to 2.54 mm and the vertical displacement at point B converges to −0.85 mm. The displacement contours at steady state are illustrated in Fig. 7a and c. For comparison, Plate-1 is also numerically studied with the commercial software ABAQUS using the linear elastic material model and the computed results are shown in Fig. 7b and d. We can find that the MPPD solutions for the displacements agree with that obtained by ABAQUS. The strain in Plate-2 is studied with the formula expressed in Eq. (28). Since the ratio of length to diameter of Plate-2 reaches 20, Plate-2 can be considered as the classical benchmark of an infinite plate with a centerized circular hole. Thus, the strain in Plate-2 can be analytically solved. The theoretical solutions are used to verify the results obtained by the proposed plane stress model. For illustration, the strain components along the positive x axis are carefully examined. The analytical solutions for the two normal strain components εx (namely, ε11 ) and ε y (namely, ε22 ) along the positive x are written as [31]: [ ] p a2 a4 εx = 2 − (5 + ν) 2 + 3(1 + ν) 4 , (90) 2E x x

14

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

Fig. 5. Square plate with a hole.

Fig. 6. The velocity and displacement at different computational steps: (a) at point A, (b) at point B.

[ ] p a2 a4 εy = −2ν + (1 + 5ν) 2 − 3(1 + ν) 4 . 2E x x

(91)

Plate-2 is uniformly discretized and numerically investigated by the same approach applied in the last numerical test. The δ-convergence (m = 3.0, δ = 3.0 mm, 5.0 mm, 7.5 mm) and m-convergence (δ = 5.0 mm, m = 3.0, 4.0, 5.0) analysis of the numerical results is illustrated in Fig. 8 for εx and Fig. 9 for ε y . As shown in Fig. 8a and b, the MPPD solutions for εx are in good agreements with the analytical solution at x/r > 5. In Fig. 8c, the numerical results of εx converges to the analytical solution as the horizon size δ increases. In Fig. 8d, m-convergence is observed near the centered hole (1 < x/r < 1.5), but the strain computed with coarse grids at 1.5 < x/r < 5 gets closer to the analytical solution. The numerical results of ε y at points away from the centered hole are slightly higher than the analytical solution, see Fig. 9a and b. Good δ-convergence is achieved in Fig. 9a and c. Although the m-convergence is not satisfied in Fig. 9b and d, the numerical results are very close.

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

15

Fig. 7. The displacement contours of u x : (a) the MPPD solution, (b) the FEM solution, and u y : (c) the MPPD solution, (d) the FEM solution.

5.2. Double cantilever beam The applicability of the proposed MPPD model for modeling deformation in solids has been demonstrated in Section 5.1. To present its ability to investigate fracturing in solids, the classical double cantilever beam (DCB) test is numerically studied. The geometry of the DCB specimen is shown in Fig. 10. The dimensions of the specimen are 0.2 m in length and 0.02 m in height. The length of the pre-crack is a0 = 0.08 m. The material properties are taken as: ρ = 7850 kg/m3 , E = 200 GPa and G 0 = 1000 J/m2 . The Poisson ratio is assumed to be ν = 0.3. The plane stress model in Section 3.3 is applied. The failure criterion given in Section 4 is used to analyze the failure process. The critical micro-potential in Eq. (89) is determined as wc = 3G 0 /(2δ 3 t), and one can refer to [32] for more details regarding on the critical value. Cracking is considered to occur once the local damage at some material point exceeds 0.4. In the numerical solution, the specimen is discretized by uniformly distributed material points with a grid spacing of 0.3 mm. The horizon size is set to be δ = 1.2 mm. The pre-crack is introduced by removing the bonds crossing the pre-crack. Displacement controlled loading u is applied at two loading points A and B. The incremental displacement at each step is set to be ∆u = 1 × 10−6 mm. The explicit time integration scheme is adopted. The time step is ∆t = 1 × 10−7 s and the total computation time is 0.2 s. At each step, we calculate the micro-potential in deformed bonds directly with Eq. (62) and use the critical micro-potential to judge whether bonds break or not. If one bond breaks, the stored potential is released. Otherwise, the potential will be taken to compute the NESED through the discretized form of Eq. (32). Besides, the kinetic energy of each discretized point is computed. The external work at each step is determined by the product of the reaction force and the incremental displacement. The total released energy is computed by the addition of energy released at each step. Simulation results are output every 200 steps. The load–displacement curve of the DCB test is shown in Fig. 11a. Due to the explicit time integration scheme, the computed reaction force oscillates especially in the process of crack growing. The reliable Smoothness Prior Method (SPM) developed by Kitagawa and Gersch [33] is used to extract the data trend from the oscillation series. The smoothed MPPD solution is illustrated in Fig. 11a. In order to facilitate analysis, all the data analyzed hereafter is smoothed by the SPM. The p–u relations for different discretizations are illustrated in Fig. 11b. As the applied displacement increases, the reaction force increases linearly to the peak value approximately. At the post-peak stage,

16

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

Fig. 8. Convergence analysis of the strain εx along the positive x axis: (a) δ-convergence, (b) m-convergence, (c) δ-convergence near the hole, (d) m-convergence near the hole.

reaction force computed with finer grids decreases gradually whereas the reaction force computed with the coarsest grid decreases sharply to zero. The difference of the p–u curves at the stage of post-peak is due to the different failure modes. For finer grids, the crack nucleats at the pre-crack tip and propagates rectilinearly along the pre-peak direction see Fig. 12a. However, the nucleated crack for the coarsest grid propagates perpendicularly and leads to the failure of the DCB specimen quickly, see Fig. 12b. The energy analysis of the DCB test with the finest grid is illustrated in Fig. 13. The total kinetic energy is very small in the whole simulation process, which reveals that the quasi-static loading condition is well satisfied. At the stage before the onset of cracking, the strain energy dominates and the total energy agrees well with the total external work. As the crack propagates, released energy increases quickly and the external work is greater than the total energy in the DCB. The difference between the external work and the total energy may be attributed to the huge oscillation of the reaction force at the post-crack stage, see Fig. 11. In Fig. 14a, the MPPD solution for the total released energy is compared with that computed by Griffith’s theory [34]. The MPPD solution is greater than Griffith’s prediction. In the simulation, the crack tip is determined by the farmost point with the local damage greater than 0.4, and thus the calculated crack growth length with such cracking assumption tends to be smaller than the actual crack growth length. As illustrated in Fig. 14b, a small amount of energy has been released as ∆a = 0.0 mm. The difference between the crack growth lengths may result in the difference between the released energy predicted by the MPPD model and Griffith’s theory. The critical loads pc and corresponding applied displacements u c for different pre-crack lengths a0 are investigated in the following. In the simulation, we determine u c as the applied displacements at which the total released energy turns out to be positive. pc are determined by the values of the smoothed p–u curves at u c . For

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

17

Fig. 9. Convergence analysis of the strain ε y along the positive x axis: (a) δ-convergence, (b) m-convergence, (c) δ-convergence near the hole, (d) m-convergence near the hole.

Fig. 10. Double cantilever beam (DCB) specimen.

comparison, theoretical predictions of pc and u c are given as [35]: ( )−3 ( ) ( a )−1 Pc uc h a0 2 0 =ι , = 32ι , (92) E Bh L h L L ( )1/2 G0 h in which ι = 96Eh . L As depicted in Fig. 15, although the computed displacements are slightly greater than the theoretical values, the MPPD solutions for critical loads agree well with the theoretical predictions. 6. Conclusions A two-dimensional micro-potential based Peridynamic formulation capable of explicitly revealing the underlying microscopic mechanism is proposed in the OSPD framework to investigate deformation and fracturing in solids.

18

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

Fig. 11. The load–displacement curves of the numerical DCB test obtained with (a) dx = 0.3 mm and (b) three different discretizations: dx = 0.3 mm, dx = 0.4 mm and dx = 0.5 mm.

Fig. 12. The crack growth paths predicted with (a) finer grids (dx = 0.3 mm, dx = 0.4 mm) and (b) the coarsest grid (dx = 0.5 mm).

Fig. 13. Energy analysis of the numerical DCB test.

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

19

Fig. 14. The variation of the released energy for different crack growth lengths.

Fig. 15. The critical loads and displacements for different lengths of the pre-crack.

Square measure of the bond length change is applied. A formulation of PD strain components for plane stress and strain problems is developed to decompose the bond length change into two parts, one resulting from the volumetric PD strain and the other from the deviatoric PD strain. At the microscale, the micro-potential generated in a deformed bond depends on the two bond length change components. At the macroscale, the nonlocal elastic strain energy density at a material point is computed through the integral of the micro-bond potential over the neighboring horizon. A general constitutive relation is derived by the Frechet derivative of the nonlocal elastic strain energy density. A failure criterion is developed based on the totally stored micro-potential. To demonstrate the reliability of the general model, several specific micro-potential functions and their corresponding constitutive relations are discussed. Numerical examples further illustrate the ability of the proposed model to model deformation and fracturing in solid materials. Acknowledgments This research was partially supported by the National Science Foundation (Grant No. 51804203) and this support is gratefully appreciated. The first author would like to thank Yong Zhang of Shanghai Jiao Tong University for his kindful help in the numerical simulation and Dr. Fei Han of Dalian University of Technology for useful discussions.

20

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

Appendix A. Equivalence of the generalized nonlocal strain components with the classical ones For homogeneous deformation, the PD strain components in Eq. (28) can be computed with Eqs. (24) and (25) as: εiPj D =

1 A

∫

ω⟨ξ ⟩

H

∆ ξ2

( 2

) ξi ξ j 1 − δi j d Sξ = L i jkl E kl − δi j Nkl E kl , ξ ·ξ 2

(A.1)

where ∫ ∫ 1 4 L i jkl = ω⟨ξ ⟩n i n j n k n l d Sξ , Nkl = ω⟨ξ ⟩n k n l d Sξ , (A.2) A H A H in which i, j, k, l = {1, 2}, n 1 = cos(θ ) and n 2 = sin(θ). In the radial coordinate system, Eq. (A.2) can be determined as: ∫ ∫ 2π 2 2π 1 L i jkl = n i n j n k n l dθ, Nkl = n k n l dθ, (A.3) π 0 2π 0 and it is easy to prove: L 1111 = L 2222 = 3/2, L 1122 = L 1212 = L 1221 = L 2112 = L 2121 = L 2211 = 1/2 and N11 = N22 = 1/2. Thus, the PD strain components in Eq. (A.1) are: PD ε11 PD ε22 PD ε12 PD ε21

= = = =

L 1111 E 11 + L 1122 E 22 − (N11 E 11 + N22 E 22 ) = E 11 , L 2211 E 11 + L 2222 E 22 − (N11 E 11 + N22 E 22 ) = E 22 , L 1212 E 12 + L 1221 E 21 = E 12 , L 2112 E 12 + L 2121 E 21 = E 21 .

(A.4) (A.5) (A.6) (A.7)

Appendix B. Calibration of parameters in the specific constitutive relation-3 Equivalent strain energy density technique is applied to establish the relationship between the material constants in the specific constitutive relation and those in the classical continuum mechanics. By substituting Eq. (67) into Eq. (32), the NESED is computed as: ∫ { [ ]} 1 − (1 + ω1 ⟨ξ ⟩|∆v |) exp −ω1 ⟨ξ ⟩|∆v | − (ω2 ⟨ξ ⟩∆d )2 d Vξ , (B.1) W PD = φ H

in which ω1 ⟨ξ ⟩ = γ ε1ξ 2 and ω2 ⟨ξ ⟩ = ζ ε1ξ 2 . c c As the deformation that the isotropic elastic solid undergoes is infinitesimal: [ ] exp −ω1 ⟨ξ ⟩|∆v | − (ω2 ⟨ξ ⟩∆d )2 ≈ 1 − ω1 ⟨ξ ⟩|∆v | − (ω2 ⟨ξ ⟩∆d )2 . Inserting Eq. (B.2) into Eq. (B.1) leads to: ∫ [ ] (ω1 ⟨ξ ⟩|∆v |)2 + (ω2 ⟨ξ ⟩∆d )2 d Sξ . W P D = φt

(B.2)

(B.3)

H

Assuming the solid ⎡ ε 0 [E] = ⎣ 0 ε 0 0

undergoes homogeneously isotropic deformation, the Green strain tensor is written as: ⎤ 0 0 ⎦, (B.4) ε33

2ν ε for plane stress problems and ε33 = 0 for plane strain problems. in which ε33 = − 1−ν The bond length change components resulting from the strain in Eq. (B.4) are computed as: { 4(1−2ν) { 2(1+ν) εξ 2 plane stress εξ 2 plane stress 3(1−ν) 3(1−ν) ∆v = ∆ = d 4 2 εξ 2 plane strain, εξ 2 plane strain. 3 3

As a result, the NESED in Eq. (B.3) can be calculated as: { 2 16C(1−2ν)2 2 ε + 4D(1+ν) ε2 plane stress PD 9(1−ν)2 9(1−ν)2 W = 16C 2 ε + 4D ε2 plane strain, 9 9

(B.5)

(B.6)

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

21

in which C =φ

∫ H

ω12 ⟨ξ ⟩ξ 4 d Vξ ,

D=φ

∫ H

ω22 ⟨ξ ⟩ξ 4 d Vξ .

The LESED under such isotropic deformation is expressed as: { 2µ(1+ν) ε2 plane stress 1−ν classical W = 2µ 2 ε plane strain. 1−2ν By equating Eq. (B.6) with Eq. (B.8), we can obtain: { 8C(1 − 2ν)2 + 2D(1 + ν)2 = 9µ(1 − ν 2 ) plane stress 9µ plane strain. 8C + 2D = 1−2ν

(B.7)

(B.8)

(B.9)

As the solid material deforms without any volumetric change, i.e. ε11 + ε22 = 0, the bond length change components are calculated as: ∆v = 0,

∆d = 2ξ T Eξ = 2εi j ξi ξ j .

(B.10)

Inserting Eq. (B.10) into Eq. (B.3) yields: W P D = DL i jkl εi j εkl ,

(B.11)

in which L i jkl is given in Eq. (A.3). Since ε11 + ε22 = 0, Eq. (B.11) can be rewritten as: W P D = Dεi j εi j .

(B.12)

In such case, the LESED is determined as: W classical = µεi j εi j .

(B.13)

Equating Eq. (B.11) with Eq. (B.13) leads to: D = µ. Inserting Eq. (B.14) into Eq. (B.9) yields: { µ(1+ν)(7−11ν) plane stress 8(1−2ν)2 C= µ(7+4ν) plane strain. 8(1−2ν)

(B.14)

(B.15)

References [1] S.A. Silling, Reformulation of elasticity theory for discontinuities and long-range forces, J. Mech. Phys. Solids 48 (1) (2000) 175–209. [2] S.A. Silling, M. Epton, O. Weckner, J. Xu, E. Askari, Peridynamic states and constitutive modeling, J. Elasticity 88 (2) (2007) 151–184. [3] S.A. Silling, R. Lehoucq, Peridynamic theory of solid mechanics, in: Advances in Applied Mechanics, Vol. 44, Elsevier, 2010, pp. 73–168. [4] S.A. Silling, R.B. Lehoucq, Convergence of peridynamics to classical elasticity theory, J. Elasticity 93 (1) (2008) 13–37. [5] Q. Du, M. Gunzburger, R. Lehoucq, K. Zhou, Analysis of the volume-constrained peridynamic Navier equation of linear elasticity, J. Elasticity 113 (2) (2013) 193–217. [6] P. Seleson, M.L. Parks, M. Gunzburger, R.B. Lehoucq, Peridynamics as an upscaling of molecular dynamics, Multiscale Model. Simul. 8 (1) (2009) 204–227. [7] E. Oterkus, E. Madenci, Peridynamic analysis of fiber-reinforced composite materials, J. Mech. Mater. Struct. 7 (1) (2012) 45–84. [8] C. Diyaroglu, E. Oterkus, E. Madenci, T. Rabczuk, A. Siddiq, Peridynamic modeling of composite laminates under explosive loading, Compos. Struct. 144 (2016) 14–23. [9] H. Zhang, P. Qiao, An extended state-based peridynamic model for damage growth prediction of bimaterial structures under thermomechanical loading, Eng. Fract. Mech. 189 (2018) 81–97. [10] S. Sun, V. Sundararaghavan, A peridynamic implementation of crystal plasticity, Int. J. Solids Struct. 51 (19–20) (2014) 3350–3360. [11] J. Luo, A. Ramazani, V. Sundararaghavan, Simulation of micro–scale shear bands using peridynamics with an adaptive dynamic relaxation method, Int. J. Solids Struct. (2018) 36–48. [12] B. Kilic, E. Madenci, Prediction of crack paths in a quenched glass plate by using peridynamic theory, Int. J. Fract. 156 (2) (2009) 165–177. [13] Y.D. Ha, F. Bobaru, Studies of dynamic crack propagation and crack branching with peridynamics, Int. J. Fract. 162 (1–2) (2010) 229–244.

22

J. Fan, R. Liu, S. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112751

[14] Y.D. Ha, J. Lee, J.-W. Hong, Fracturing patterns of rock-like materials in compression captured with peridynamics, Eng. Fract. Mech. 144 (2015) 176–193. [15] X. Lai, B. Ren, H. Fan, S. Li, C. Wu, R.A. Regueiro, L. Liu, Peridynamics simulations of geomaterial fragmentation by impulse loads, Int. J. Numer. Anal. Methods Geomech. 39 (12) (2015) 1304–1330. [16] X. Gu, Q. Zhang, D. Huang, Y. Yv, Wave dispersion analysis and simulation method for concrete SHPB test in peridynamics, Eng. Fract. Mech. 160 (2016) 124–137. [17] H. Fan, S. Li, A peridynamics-SPH modeling and simulation of blast fragmentation of soil under buried explosive loads, Comput. Methods Appl. Mech. Engrg. 318 (2017) 349–381. [18] W. Gerstle, N. Sau, S. Silling, Peridynamic Modeling of Plain and Reinforced Concrete Structures, IASMiRT, 2005. [19] Q. Zhu, T. Ni, Peridynamic formulations enriched with bond rotation effects, Internat. J. Engrg. Sci 121 (2017) 118–129. [20] Y. Wang, X. Zhou, Y. Wang, Y. Shou, A 3-D conjugated bond-pair-based peridynamic formulation for initiation and propagation of cracks in brittle solids, Int. J. Solids Struct. 134 (2017) 89–115. [21] Y. Hu, E. Madenci, Bond-based peridynamic modeling of composite laminates with arbitrary fiber orientation and stacking sequence, Compos. Struct. 153 (2016) 139–175. [22] H. Ren, X. Zhuang, T. Rabczuk, A new peridynamic formulation with shear deformation for elastic solid, J. Micromech. Mol. Phys. 1 (02) (2016) 1650009. [23] S.A. Silling, E. Askari, A meshfree method based on the peridynamic model of solid mechanics, Comput. Struct. 83 (17–18) (2005) 1526–1535. [24] Q. Le, W. Chan, J. Schwartz, A two-dimensional ordinary, state-based peridynamic model for linearly elastic solids, Internat. J. Numer. Methods Engrg. 98 (8) (2014) 547–561. [25] Y. Zhang, P. Qiao, An axisymmetric ordinary state-based peridynamic model for linear elastic solids, Comput. Methods Appl. Mech. Engrg. 341 (2018) 517–550. [26] H. Ren, X. Zhuang, Y. Cai, T. Rabczuk, Dual-horizon peridynamics, Internat. J. Numer. Methods Engrg. 108 (12) (2016) 1451–1476. [27] H. Ren, X. Zhuang, T. Rabczuk, Dual-horizon peridynamics: A stable solution to varying horizons, Comput. Methods Appl. Mech. Engrg. 318 (2017) 762–782. [28] E. Madenci, Peridynamic integrals for strain invariants of homogeneous deformation, ZAMM Z. Angew. Math. Mech. 97 (10) (2017) 1236–1251. [29] B. Kilic, E. Madenci, An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory, Theor. Appl. Fract. Mech. 53 (3) (2010) 194–204. [30] T. Rabczuk, H. Ren, X. Zhuang, A nonlocal operator method for partial differential equations with application to electromagnetic waveguide problem, Comput. Mater. Continua 59 (2019) Nr. 1. [31] N. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity, Springer, 1977. [32] D. Dipasquale, G. Sarego, M. Zaccariotto, U. Galvanetto, A discussion on failure criteria for ordinary state-based peridynamics, Eng. Fract. Mech. 186 (2017) 378–398. [33] G. Kitagawa, W. Gersch, Smoothness Priors Analysis of Time Series, Vol. 116, Springer Science & Business Media, 1996. [34] T.L. Anderson, Fracture Mechanics: Fundamentals and Applications, CRC Press, 2005. [35] A. Abbaszadeh Bidokhti, A.R. Shahani, M.R. Amini Fasakhodi, Displacement-controlled crack growth in double cantilever beam specimen: a comparative study of different models, Proc. Inst. Mech. Eng. C 231 (15) (2017) 2835–2847.