A state-based peridynamic model for quantitative elastic and fracture analysis of orthotropic materials

Accepted Manuscript A state-based peridynamic model for quantitative elastic and fracture analysis of orthotropic materials Heng Zhang, Pizhong Qiao PII: DOI: Reference:

S0013-7944(18)30718-5 https://doi.org/10.1016/j.engfracmech.2018.10.003 EFM 6177

To appear in:

Engineering Fracture Mechanics

Received Date: Revised Date: Accepted Date:

18 July 2018 1 October 2018 3 October 2018

Please cite this article as: Zhang, H., Qiao, P., A state-based peridynamic model for quantitative elastic and fracture analysis of orthotropic materials, Engineering Fracture Mechanics (2018), doi: https://doi.org/10.1016/ j.engfracmech.2018.10.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A state-based peridynamic model for quantitative elastic and fracture analysis of orthotropic materials Heng Zhanga, Pizhong Qiaoa,b, a

State Key Laboratory of Ocean Engineering, Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, PR China b

Department of Civil and Environmental Engineering, Washington State University, Sloan Hall 117, Pullman, WA 99164-2910, USA

Abstract A new state-based peridynamic model is proposed for quantitative elastic and fracture analysis of the orthotropic materials. In this model, four independent material parameters of two-dimensional (2D) orthotropic materials in the conventional composite mechanics model are transformed into the same number of peridynamic constants that lead to no constraints on material properties, and three types of critical stretches are uniquely defined for quantitative fracture analysis of orthotropic materials. Quantitative simulations of an orthotropic lamina under stress load and the compact tension (CT) and single edge-notched tension (SENT) tests are performed to verify the respective elastic and fracture behaviors of orthotropic materials by the proposed model. The strain of lamina and the typical fracture parameters (i.e., compliance, critical load, and released energy density) are calculated and compared with available analytical and experimental data. The present peridynamic model is capable of capturing fracture characteristics (i.e., crack path, energy concentration, displacement-load



Corresponding author. E-mail addresses: [email protected] and [email protected] (P. Qiao). 1

relationship, and energy balance during crack propagation) of orthotropic materials, as demonstrated in both the CT and SENT tests. Keywords: composite mechanics; peridynamics; orthotropic materials; state-based model; fracture analysis. 1. Introduction The anisotropic materials, such as woods, ceramics and composites, are widely used for engineering structures. Prediction of the fracture behavior of the anisotropic materials and corresponding fracture analysis are significant in engineering applications. In order to understand fracture characteristics of the anisotropic materials, available numerical methods established for the isotropic materials have been extended and applied for failure analysis of anisotropic materials, such as finite element method (FEM) [1], extended finite element method [2], boundary element method [3], etc. However, these methods require quite complex laws, and they are all based on classical continuum mechanics theory, which is naturally not fit in the discontinuous fields. The peridynamic theory was first proposed by Silling [4] to overcome the discontinuous problems which involves fracture and crack propagation. Instead of differentiation, integration is employed by peridynamics in its equations of motion; thus, it can avoid the derivatives of displacement in its discontinuous fields. Furthermore, “bonds” are used between material points to describe the nonlocal interactions, and damage in structures appears naturally with the developing broken bonds; thus, no additional fracture criteria are needed as in FEM when applied to crack growth problems. Generally, the peridynamic model can be classified into “bond-based” [4] and “state-based” [5] peridynamics, and the original

2

bond-based model was first proposed but with the constraint on the material properties. The state-based model is a more general framework that employs explicit dependence of volumetric and distortional deformations for material response and thus reduces the limitations on material properties of original model. With the great capability of dealing with discontinuous situation, the peridynamic model has been successfully applied for failure predictions of isotropic materials [6–8]. Peridynamics has also been extended for fracture analysis of anisotropic materials [9–14]. Askari et al. [9] investigated the failure modes of laminated composites with a mixed constitutive peridynamic model for different bonds of fibers and matrix. Gerstle et al. [10] studied the peridynamic model of plain and reinforced concrete structures. Xu et al. [11] performed damage and failure analysis of composite laminates under biaxial loads. Oterkus et al. [12] introduced the prediction of failure in a stiffened composite curved panel with the combined finite element and peridynamics. Hu et al. [13] proposed a homogenized peridynamic composite model for dynamic fracture analysis in unidirectional fiber-reinforced composites. Ghajari et al. [14] performed the analysis of dynamic crack propagation in orthotropic media by using a continuous function of bond parameters. However, all the above models were based on the bond-based peridynamic theory, and the number of the material parameters are reduced into less number of peridynamic constants, leading to the limitations on the material properties. Madenci et al. [15] proposed a more general model for composite lamina to remove those limitations, but it was still in the discrete form and with the special form of influence function. In this paper, a new state-based peridynamic model for two-dimensional (2D)

3

orthotropic materials is proposed, and an orthotropic lamina under stress load and the compact tension (CT) and single edge-notched tension (SENT) tests are quantitatively studied for validations of the respective elastic and fracture behaviors of the new model. First, a state-based peridynamic model for 2D orthotropic materials is proposed, three types of critical stretches are defined for fracture analysis. Then, the elastic behavior of an orthotropic lamina under uniaxial tensile load is studied with its associated convergence study, followed by the elastic analysis of the orthotropic lamina with different load conditions; while the fracture behaviors of the CT and SENT tests are analyzed. Crack predictions of the CT specimens are performed with different pre-crack lengths or anisotropy (or fiber) orientations. The different values of critical stretches s2 and s1 are discussed to investigate their effect on fracture analysis.

2. Orthotropic model In this section, a theoretical background of state-based peridynamic model for isotropic material is first introduced, followed by the proposed state-based peridynamic model for two-dimensional (2D) orthotropic materials. 2.1. A brief review of peridynamic model for isotropic materials The equation of motion of the material point x in the general state-based peridynamic model can be expressed as [16]:

  x  u  x, t   

x

  x, t 

x  x    x, t  x  x  dVx  b  x, t 

(1)

where ρ is the mass density, u is the displacement of point x at time t, x  is another material point in the neighborhood Hx of x, and b  x,t  is the body force density of point x. As shown in Fig. 1,   x  x is the bond vector, and   x,t  and   x,t  are the force 4

vector states of points x and x  , respectively.

  x,t  x  x

x

  x ,t  x   x



x δ Hx

Fig. 1. Ordinary state-based peridynamic model.

In the ordinary state-based peridynamic model, the general strain energy density for the linear elastic isotropic material takes the form of [5]:

k 2  d d W  e e 2 2



where 

is the volume dilatation, e

d



(2)

is the deviatoric extension state,  is the influence

function, and k and α are the peridynamic constants, respectively. And the dot product



is

defined in [5]. For the three-dimensional (3D) and two-dimensional (2D) cases, those parameters are expressed as [17]:

x  e 15  d   3 q , e  e   x 3, k  k ,   q ;      2  x  e , ed  e   x 2, k  k ,   8 ;  q q

3D (3)

2D

And the scalar force state t for the linear elastic materials takes the forms of [17]:

   x 15 d 3k  q  q e  t  2k   x  8 e d  q q

3D (4)

2D

where x is the bond length scalar state of bond  , e is the extension scalar state, q is the

5

weighted volume defined as  x  x , μ is the shear modulus, and k  is the bulk modulus that can be expressed in terms of the Young’s modulus E and the Poisson’s ratio v as:

 E 3D  3 1  2 v     E k   Plane stress 2 1  v     E Plane strain   2 1  v 1  2v 

(5)

2.2. Peridynamic model for orthotropic materials 2.2.1 Peridynamic strain energy density The form of the peridynamic energy density for 2D orthotropic materials can be expressed as:

k 2  1 1 d d W   e  e  c1 1 e   e  c2  2 e   e 2 2 2 2





where k, α, c1 and c2 are the peridynamic constants to be found. θ and e

(6) d

are the volume

dilatation and deviatoric extension state, respectively, defined as:

 2 where the dot product



x  e q

e  e  x 2 d

,

(7)

is defined in [5].  1 and  2 are the influence functions that

take the forms of:

1       ,       ,         3    2       ,        ,     2 2  



6





(8)

y x 

2

 

1 x δ

x Hx

Fig. 2. State-based peridynamic model for 2D orthotropic materials.

As shown in Fig. 2,   x  x is the vector state,  

is the scalar state with the

angle of the bond  in x-y coordinate, and φ is the orientation of the first principal axis to the x-axis. The Dirac-delta function    ,   can be expressed as:

0       1    

   ,    

(9)

So  1 is equal to the influence function ω when the bond  is parallel to the 1-principal material axis (     or       ), while  2 is equal to ω when  is parallel to the 2-axis; otherwise,  1 and  2 are equal to 0. Thus, using Eqs. (8) and (9), the dot product with those special influence functions in Eq. (6) can be calculated as:

1 e   e      ,      ,     e2 dA  L  e2 d 1

  3  2 e   e      ,       ,     e2 dA  L  e2 d 2 2  







(10)

2

where L1 and L2 are the anisotropic effect lines in domain H parallel to the 1 and 2 axes, respectively. In addition, similar to the definition of q, the weighted volume m1 and m2 are defined as m1  1 x   x and m2  2 x   x , and they can be regarded as the weighted volume of the 1- and 2-axes, respectively.

7

In the present orthotropic peridynamic model, there are three types of “bond” between material points, i.e.,  0 is the bond of equivalent isotropic model corresponding to the first two parts of strain energy density W; 1 and  2 are the additional bonds along the anisotropic effect lines L1 and L2 corresponding to the third and fourth parts of W, respectively. In summary, the total strain energy density of the 2D orthotropic materials in Eq. (6) can be regarded as the summation of the energy from the equivalent isotropic model part and the respective anisotropic effect energies in its 1- and 2-principal material axes. 2.2.2 Peridynamic constants Since the strain energy densities of point in peridynamics and classical mechanics are equal, the unknown peridynamic constants in Eq. (6) can be obtained by equalizing two energies. From classical orthotropic material theory, the relationship of stress vs. strain in the 1-2 plane can be described as:  11  C11     22   C12    0  12  

C12 C22 0

0   1    0    2   C66    12 

(11)

where Cij is the stiffness matrix. The stress-strain relationship in arbitrary orientation is:

 x   C11 C12     y   C12 C 22     xy  C16 C 26

C16    x    C 26    y    C 66   xy 

where C ij is transformed from Cij with anisotropy (fiber) orientation φ as:

8

(12)

C11  C11 cos 4   2  C12  2C66  sin 2  cos 2   C22 sin 4 

C12  C12  sin 4   cos 4     C11  C22  4C66  sin 2  cos 2  C 22  C11 sin 4   2  C12  2C66  sin 2  cos 2   C22 cos 4 

(13)

C16   C11  C12  2C66  sin  cos3    C12  C22  2C66  sin 3  cos  C 26   C11  C12  2C66  sin 3  cos    C12  C22  2C66  sin  cos3  C 66   C11  C22  2C12  2C66  sin 2  cos 2   C66  sin 4   cos 4   The strain energy density of a material point under the strain state  i is:

1 Wc  C ij  i  j 2 where i and j = 1,2,6 and  6   12 .  1 ,  2 ,  12 

T

(14) T

can also be expressed as  x ,  y ,  xy  . As

shown in Fig. 2, the extension scalar state e of bond  can be expressed with the strain  x ,  y ,  xy 

T

and the direction scalar state  

as:

e   cos2  x  sin 2  y  sin  cos xy  x

(15)

As shown in Fig. 3, three kinds of strain states are considered to obtain the unknown peridynamic constants. Using Eqs. (7) and (15), the respective values of e, θ and e

d

for the

three typical strain states are reported in Table 1. L ζ

y

2

2

y

2

1

1

W

W

x

L

(a)

(W/2)ζ

1

W

W+Wζ

x

x L

L

L

y

L+Lζ

L+Lζ

(b)

(c)

Fig. 3. Three types of strain states: (a) simple shear, (b) x-directional uniaxial stretch, and (c) biaxial stretch.

9

(1) Simple shear: By substituting the strain state of simple shear in Table 1 into Eq. (14), the classical strain energy density of a material point under the simple shear strain state is obtained as:

1 Wc  C 66 2 2 By substituting the values of e, θ and e

d

(16)

for the simple shear strain case in Table 1 into

Eq. (6) and using the influence function of Eq. (8), the peridynamic energy density under the simple shear is expressed as:

WPD 

1 1 1  q 2  c1m1 sin 2  cos2  2  c2 m2 sin 2  cos2  2 16 2 2

(17)

where m1 and m2 are the weighted volumes in the 1- and 2-principal material axes, respectively. By equalizing the classical and peridynamic energy densities in Eqs. (16) and (17) and using the definition of C 66 in Eq. (13), the following equations are obtained as:

 q 8  C66   c1m1  c2 m2  C11  C22  2C12  4C66

(18)

(2) Uniaxial stretch in x-direction: The strain energy density of a material point under uniaxial stretch is:

1 Wc  C11 2 2 By using the values of e, θ and e

d

(19)

for uniaxial stretch strain case in Table 1, the

peridynamic energy density under uniaxial stretch is given as:

1 1 1 1 WPD  k 2   q 2  c1m1 cos 4  2  c2 m2 sin 4  2 2 16 2 2

(20)

By equalizing the classical and peridynamic energy densities in Eqs. (19) and (20) and using the definition of C11 in Eq. (13), the following relationships are obtained as: 10

k   q 8  c2 m2  C22   2c2 m2  2C12  C22  4C66  c m  c m  C  C  2C  4C 11 22 12 66 1 1 2 2

(21)

(3) Biaxial stretch: The strain energy density of a material point under biaxial stretch is:

Wc 





1 C11  2C12  C 22  2 2

Again, by using the values of e, θ and e

d

(22)

for biaxial stretch strain case in Table 1, the

peridynamic energy density under biaxial stretch is given as:

1 1 WPD  2k 2  c1m1 2  c2 m2 2 2 2

(23)

By equalizing the classical and peridynamic energy densities in Eqs. (22) and (23) and using the definitions in Eq. (13), the following relation is obtained as:

4k  c1m1  c2 m2  C11  2C12  C22

(24)

In summary, though there are six equations for four unknown peridynamic parameters in Eqs. (18), (21) and (24), only four equations are independent. Thus, the values of peridynamic parameters are obtained uniquely in term of lamina stiffness components as:

k  C12  C66   8C    66 q    c  C11  C12  2C66  1 m1  C22  C12  2C66  c2  m2 

(25)

Specially, according to the expressions in Eq. (25), when the orthotropic material parameters satisfy: C22 < C12 + 2C66 < C11, the peridynamic constants meet the conditions of c1 > 0 and c2 < 0. According to the peridynamic strain energy density expression in Eq. (6),

11

when c2 < 0, the fourth parts of W is negative, which means that the additional bonds  2 along the anisotropic effect lines L2 have the negative effect on its strain energy density. However, it is acceptable when the point have the intact horizon. Even though the fourth part of the W is negative, the whole expression W is positive for any given deformation state.

Table 1. Three types of strain states Cases

Strain states

E

θ

ed

Simple shear

 xy   ,  x   y  0

sin cos x

0

sin cos x

Uniaxial stretch

 x   ,  y   xy  0

cos2  x



1 cos 2 x 2

Biaxial stretch

 x   y   ,  xy  0

x

2

0

2.2.3 Peridynamic force state In the ordinary model [5], the scalar force state t can be calculated by taking the Frechet derivative of peridynamic energy density W with respect to e, as teW . Taking the Frechet derivative of function of Eq. (6), the scalar force state t for the 2D orthotropic materials is thus obtained as:

t

2 x d k   e  c11 e  c22 e q

(26)

where k, ??, c1 and c2 are given in Eq. (25).  and e are defined in Eq. (7). d

For isotropic materials, the stiffness matrices in the cases of plane stress and plane strain can be expressed as:

12

    1 v 0    E  v 1 0   2  1 v   1 v  0 0    2  Cij     1  v v   E  v 1 v  1  v 1  2v     0  0   

Plane stress

 0   0  Plane strain 1  2v   2 

(27)

Substituting Eq. (27) into Eq. (25), the values of peridynamic constant for isotropic materials can be rewritten as:

 k  k    8 q   c 0  1  c2  0

(28)

By using the values in Eq. (28), the strain energy density in Eq. (6) and scalar force state in (26) of the present model can be reduced to those of isotropic materials, as given in Eqs. (2) and (4), respectively. In summary, the general forms of the strain energy density and scalar force state of state-based peridynamics for 2D orthotropic materials are given in Eqs. (6) and (26), respectively. Unlike the bond-based orthotropic model [13,14] where four independent orthotropic material parameters were reduced to fewer peridynamic constants, thus leading to the existing restrictions on material properties, the relationships of four peridynamic constants with four orthotropic material parameters are for the first time established in the present new model (see Eq. (25)), which thus captures the general characteristics of orthotropic materials. Also, unlike the peridynamic model for composite lamina given in [15], which was in the discrete form and with the special influence function, the present new model 13

developed is a general state-based one for the orthotropic materials and suitable for any reasonable influence function. In addition, as a validation, the proposed peridynamic model for orthotropic materials is simplified and reduced to the peridynamic model for isotropic materials with the typical material parameters of isotropic cases. 2.3. Peridynamic fracture analysis model for orthotropic materials The suitable peridynamic damage model is necessary and important for analysis of fracture problems. Silling and Askari [18] first proposed the critical stretch damage model, which was based on the bond-based peridynamics. Most recently, the authors [17] proposed a state-based peridynamic model for quantitative fracture analysis of isotropic materials and for the first time provided the general relationship of critical stretch for bond failure with respect to the critical energy release rate of the isotropic materials. In this study, the fracture model for isotropic materials developed by the authors [17] is extended to deal with the case of the 2D orthotropic materials. Similar to the one presented in [16], the strain energy density of peridynamic model for orthotropic material in Eq. (6) can be rewritten as:

1  1 1 W  a 2   e   e  c1 1 e   e  c2  2 e   e 2 2 2 2

(29)

where a is defined as:

a  k 

q 4

(30)

The bond stretch scalar state s is defined as: s = (e- x)/x [18]. Since there are three types of bonds with typical bond parameters for orthotropic materials, three types of critical stretches are defined: s0 is the critical stretch of bond  0 for equivalent isotropic model part,

14

and s1 and s2 are the critical bond stretches of additional bonds 1 and  2 in the 1- and 2-principal material axes, respectively. When the bond stretch of  i (i = 0, 1, 2) grows beyond its respective critical stretch si (i = 0, 1, 2), the corresponding typical bond is broken, and the respective history-dependent scaler μi (i = 0, 1, 2) changes from 1 to 0 as:

1 if s  t , i   si for all t   t , i  0, 1, 2 otherwise 0

 i i   

(31)

As shown in Fig. 4, as the crack grows along the angle of φ (anisotropy or fiber orientation), the bond stretch of  i across the crack all arrives at the critical stretch si and the extension scalar e  si x . The released energy of a material point is equal to its loss of the strain energy during crack propagation. For the equivalent isotropic model part, the bond stretch of the broken  0 is equal to s0 and the volume dilatation caused by all those broken bonds in domain H1 is:

s  2

 xe q

 2

1



 x 2 s0 dV



 x 2 dV

 2 ( z ) s0

(32)

where λ (z) is the ratio of the weighted volume from H1 to H domain, defined as:

  ( z)  

1 

 x 2 dV  x 2 dV

 

1

 x 2 dV (33)

q

y x

2 H1

1



x z

δ o

Fig. 4. The total energy absorbed of point x in the domain H1 as crack propagates along the

15

anisotropy angle of φ.

As the crack grows, the additional bonds 1 and  2 across the crack is equal to s1 and s2, respectively, and the released energy of the orthotropic material thus varies with the crack orientation because of those additional bonds along the principal axes. By substituting Eq. (32) into the strain energy density in Eq. (29) and using Eq. (33), the released energy density of a point because of the crack propagation takes the form of:

W

a 2  1 1  s    x 2 s02 dV  c1  1 x 2 s12 dV  c2   2 x 2 s22 dV 2 2 1 2 1 2 1 1 1 1   4a ( z ) 2   q ( z )  s0 2  c1m1  1 ( , z ) s12  c2 m2  2 ( , z ) s2 2 2 2 2

(34)

where 1 ( , z ) and  2 ( , z ) are defined along the principal axes, as the ratios of the weighted volume from parts of principal axes of L1 and L 2 in H1 domain to the whole principle axes of L1 and L2, respectively,

  x dA    x d    x d ,  ( , z )  m   x dA   x d   x dA    x d    x d  ( , z )  m   x dA   x d 2

1

1

1



1

2

2

1

1

L1

2

2

L2

2



L1

2

2

2

2

L1

2

(35)

L2

2

2

L2

2

Typically, when the crack grows along the plane normal to the 1-axis (φ = 0°) that leads to

1 ( ,z ) 1 z( )and  2 ( , z )  0 ; while the crack path normal to the 2-axis (φ = 90°) leads to 1 ( , z )  0 and  2 ( , z )  2 ( z ) , where λ1(z) and λ1(z) are the typical values of line weighted volume ratios as φ = 90° and φ = 0°, respectively. Thus, the released energy densities of a point as the crack normal to the 1- and 2-axes can be, respectively, written as: 1 1  W 1   4a ( z ) 2   q ( z )  s0 2  c1m11 ( z ) s12   2 2  W 2  1  4a ( z ) 2   q ( z )  s 2  1 c m  ( z ) s 2 0 2 2 2 2   2 2

16

(36)

Thus, the total released energies at unit area normal to the 1- and 2-axes are, respectively, written as:

 dU e1  2  dW 1 dz   4a    q  s 2  c m  s 2 0 1 1 1 1 0    dU e 2  2 dW 2 dz   4a    q  s0 2  c2 m2  2 s2 2 0 

(37)

where  ,   , 1 and 1 are the parameters that are affected by the influence function defined as: 







0

0

0

0

    ( z )dz,      ( z )2 dz, 1   1 ( z )dz,  2   2 ( z )dz

(38)

The total released energies at unit area in the plane normal to the 1- and 2-axes are equal to the critical energy release rates in their respective axes, and for the mode Ι type crack, they are expressed as: 2 2   GIc1  dU e1   4a    q  s0  c1m11s1  2 2  GIc 2  dU e 2   4a    q  s0  c2 m2  2 s2

(39)

Typically, for the isotropic materials, using Eq. (28), the form of Eq. (39) can be reduced to: s0  GIc

 4k   8     8 

(40)

which is equal to the form of critical stretch of isotropic material under plane stress and strain conditions [17]. Because only two equations are given for three unknown parameters in Eq. (39), si does not have the unique solution. As aforementioned (see the discussion after Eq. (25)), when c2 < 0, the fourth parts of W is negative, which means that the additional bonds  2 have the negative effect on its strain energy density. Thus, when the bonds  2 are broken at the time the stretch of  2 reach its critical stretch s2, the responding strain energy density consequently increases and the released energy decreases. It is physically unacceptable that 17

the released energy inversely decreases during crack propagation. Only when s0 = s2, the bonds  2 are broken together with amount of bonds  0 , and the released energy will increase with the broken bonds at any time. In addition, because for the unidirectional composite materials, the anisotropic effect in the 2-principal material axes depends on the matrix, it is reasonable that the critical stretches meet the condition of s0 = s2. Thus, when s0 = s2 is adopted for the two dimensional (2D) orthotropic materials, the critical bond stretch si can be uniquely obtained from Eq. (30). In addition, the point variable parameters di (i = 0, 1, 2) as the volume weight of broken bonds are defined for the bond  i , respectively, as:

d 0  x, t   1 



H

 0 dV



H

dV

, d1  x, t   1 

  dV , d x, t  1    dV   dV    dV L1

1

L2

2

L1

2

(41)

L2

where μi is the history-dependent scaler of point x defined in Eq. (31). Thus, the variable parameters di embody the broken state of bonds in orthotropic peridynamic model, and it can be used to trace crack path during simulation. Moreover, if without extra statement, the typical d0 is used to trace the crack path in this paper. In summary, the state-based peridynamic fracture analysis model for 2D orthotropic materials is developed, and three types of critical stretches are defined and calculated from the critical energy released rates in Eq. (39). As a validation, the proposed fracture analysis model for orthotropic materials can be reduced to the damage model for isotropic materials in Eq. (40) [17], when the parameters of isotropic materials are used. In addition, the released energy density of orthotropic materials is given in Eq. (35) to quantitatively trace the released energy of a point during the crack propagation, and the sum of the released energy of all

18

points is equal to the total incremental surface energy. 3. Numerical method Similar to the study in [18], the system is discretized into the finite material nodes, and each node has a finite volume. The equation of motion in Eq. (1) can thus be rewritten as:

i ui 

 {T [ x , t ] i

ji

x j  xi  T [ x j , t ] xi  x j }V j  bi

(42)

where T is defined in Eq. (26) and Vj is the volume of node j. For the 2D problems with uniform grid, V j  h x 2 , where x is the uniform grid size and h is the constant thickness of 2D plate. The summation of Eq. (42) is taken over in the horizon Hi of node i as x j  xi   , where δ is the horizon size as   mx . y

2

2γ 1

o

x

Fig. 5. The discrete peridynamic model for orthotropic materials.

As shown in Fig. 5, after the grid discretization, the nodes along the 1- and 2-principal material axes need to be found. Because the discretization grid may not align with the principal orientations, the angle tolerance  is used to find the node for the additional bonds [13]. Unlike the study in [13], where the bonds are explicitly divided into “fiber bonds” and

19

“matrix bonds” and the material node belongs to unique kind of bond, the point along the principal axes is connected with two kinds of bonds in the present model: the one corresponding to the equivalent isotropic model and the other of additional bond corresponding to the anisotropic effect. For example, as shown in Fig. 5, for the typical value of φ and the angle tolerance  , the blue and red nodes along the principal axes contribute not only to the equivalent isotropic part but also the associated anisotropic part of strain energy, and they connect to the central node with two kinds of bonds:  0 and 1 or  0 and  2 , respectively. In addition, the angle tolerance  used is 5° in this study. Even though the number of nodes for the additional bonds connected to the central node varies with the orientation of the principal axes, the strain energy density of point in Eq. (6) is independent of the principal orientation because of the characteristic of state-based peridynamic model. Thus, the scaling work of bonds for the strain energy density equivalence in [13] is not needed. In this study, two kinds of convergence studies [19] are performed: m-convergence and δ-convergence. The m-convergence is applied when the value of  is fixed and the value of m increases, while δ-convergence is evaluated when m is fixed and  decreases. However, the nodes around boundary or discontinuous crack face do not have a full horizon, which leads to the skin effect [20]. The effective material properties of boundary nodes are different from those in bulk, because the skin effect depends on the horizon size and influence function [20]. In this paper, the denser mesh size from convergence studies is adopted, and the typical conical influence function    / x is utilized to reduce the skin effect. Moreover, the adaptive refinement [21] or the surface correction [20] can be the

20

reasonable methods to further solve the skin effect. 4. Applications In this section, the elastic behavior of an orthotropic lamina under static stress loads is first analyzed to verify its accuracy in elastic analysis, followed by the simulations of both compact tension (CT) and single edge-notched tension (SENT) tests to demonstrate its application to quantitative fracture analysis of orthotropic materials. 4.1 Orthotropic lamina The geometries of the orthotropic lamina are given in Fig. 6. As shown in Fig. 6, under T

the stress of  x ,  y , xy  , the strain of  x ,  y ,  xy 

T

for various angles of anisotropy

(fiber) is obtained from the peridynamic simulation and compared to its theoretical values determined from:

  x   S 11      y    S 12  xy   S 16   

S 12 S 22 S 26

S 16   x    S 26   y   S 66   xy 

(43)

where S ij is the compliance matrix of the lamina, which is the inverse of the stiffness matrix shown in Eq. (13). For the case of plane stress condition problems, Cij = Qij are given as: Q11 

E1 E2 v E , Q22  , Q12  12 2 , Q66  G12 1  v12v21 1  v12v21 1  v12v21

(44)

where E1 and E2 are the elastic Young’s moduli in the 1- and 2-principal material directions, respectively, v12 and v21 are the Poisson’s ratios, and G12 is the shear modulus in 1-2 plane. By substituting Eq. (44) and the material values in Table 2 into Eq. (13), the stiffness matrix C ij for anisotropy is obtained, so is the compliance matrix S ij by inverting C ij . Then, the

21

strain components under any stress state are obtained from Eq. (43) as the theoretical values. In this example, the material used is the IM-7/977-3 carbon fiber composite [22], and the material properties are reported in Table 2. The plane stress condition is considered, and the constant value of lamina thickness of h = 1 mm is used. In addition, seven different anisotropy (fiber) angles φ = 0°,15°, 30°, 45°, 60°, 75° and 90° are investigated. The adaptive dynamic relaxation (ADR) method [23] is utilized to obtain the final stable results. Since the strain is not defined in the peridynamic model and the displacement distributions of composite lamina is linearly smooth under the stress load, the average strain is defined. Typically, the displacements of points (Fig. 6): A (10 mm, 20 mm), B (30 mm, 20 mm), C (20 mm, 10 mm), and D (20 mm, 30 mm) are used to calculate the average values of strain as:

 u x (B)  u x (A) ,  x  20 mm   u y (D)  u y (C) ,  y  20 mm   u y (B)  u y (A) u x (D)  u x (C)   xy  20 mm 20 mm 

(45)

where ux(A), ux(B), ux(C), and ux(D) are the x component displacements of points A, B, C, and D, respectively; while uy(A), uy(B), uy(C) and uy(D) are the respective y component displacements. In the following parts, the orthotropic lamina under uniaxial tensile load is first investigated to validate the convergence of the proposed model by the typical m-convergence and δ-convergence. The other three applied loads (i.e., the simple shear, y-directional tension, and biaxial tension) are then considered to demonstrate suitability and robustness of the new

22

model for complex load conditions.

τxy 2

1

40 mm

D

σx

A

B

τxy

C y

x

40 mm

o

σy Fig. 6. Orthotropic lamina under stress load. Table 2. Elastic properties of IM-7/977-3 carbon fiber composite [22] E1

E2

G12

(GPa)

(GPa)

(GPa)

164.3

8.977

5.02

v12

ρ (kg/m3)

0.32

1603

4.1.1 Orthotropic lamina under uniaxial tensile load As shown in Fig. 6, the plate under uniaxial tensile load along the x-direction is first considered for the convergence study, and the typical stress state takes the form of

 x  1MPa,  y   xy  0 . First, a fixed horizon size of δ = 2 mm and different values of m = 4, 6, 8 are used for m-convergence test. The strain components of orthotropic lamina under x-directional tensile load with different values of m are shown in Fig. 7. As shown in Fig.7, the peridynamic solution converges to the theoretical results as m increases. Typically, when φ = 0° and 90°,

23

the simulation model accurately predicts the strain values with even small value of m, within difference of 1%. As for the other values of φ, the strain components of εx and γxy converge to the theoretical values when m increases to 8 with the maximum difference of 6.4%. The strain component of εy is more difficult to converge to its theoretical value, especially when φ is close to 45°, with the maximum difference of 7.9%. However, the strain in the Poisson’s direction (i.e., in the direction normal to the loading) is much smaller than the other strain components, and it slightly affects the whole strain response. 12 10 8 6 4

m=4 m=6 m = 8

-1.0 -1.5 -2.0

2 0

Classical

-0.5 y (× 10-3)

x (× 10-3)

0.0

Classical m=4 m=6 m = 8

0

15

30

45 

60

75

-2.5

90

0

15

30

(a)

45 

60

75

(b)

2

Classical m = 4

xy (× 10-3)

0

m=6 m=8

-2 -4 -6 -8

0

15

30

45 

60

75

90

(c) Fig. 7. Strain components of orthotropic lamina under x-directional tensile load with δ = 2 mm and different m: (a) εx, (b) εy, and (c) γxy.

24

90

Then, the fixed value of m = 8 is selected to perform the δ-convergence test with the varying values of δ = 4 mm, 2 mm, and 1 mm, respectively. The strain of orthotropic lamina under x-directional tensile load with different values of δ are shown in Fig. 8. As shown in Fig. 8, the simulation results converge to the theoretical values as the horizon δ decreases. Typically, the peridynamic model provides accurate solutions for εx and γxy with the maximum difference of 4.2%; while the maximum difference of εy is 6.1% in the Poisson’s direction (i.e., in the direction normal to the loading). 12

0.0  Classical  = 4 mm  = 2 mm  = 1 mm

10

y (× 10-3)

x (× 10-3)

8 6 4

-1.0 -1.5 -2.0

2 0

 Classical  = 4 mm  = 2 mm  = 1 mm

-0.5

0

15

30

45

60

75

-2.5

90

0

15

30

45





(a)

(b)

60

75

90

2  Classical  = 4 mm  = 2 mm  = 1 mm

xy (× 10-3)

0 -2 -4 -6 -8

0

15

30

45 

60

75

90

(c) Fig. 8. Strain components of orthotropic lamina under x-directional tensile load with different δ and m = 8: (a) εx, (b) εy, and (c) γxy.

25

Unlike the study in [14] in which only the strain along the loading direction was T

obtained, the complete set of strain components  x ,  y ,  xy  are obtained in the present example. Furthermore, the more numbers of typical anisotropy (fiber) angles φ are considered to validate the accuracy of the present model for arbitrary orientation. Similar to the isotropic material model in [24], in which the estimated values of Poisson’s ratio have errors because of boundary effects, the strain in the Poisson’s direction (i.e., in the direction normal to the loading) is also affected by the boundary effects in the present example. 4.1.2 Orthotropic lamina under three different load conditions Three different load cases are considered to validate capability and robustness of the proposed model for complex load conditions. From the convergence study above, the values of δ = 2 mm and m = 8 are chosen to obtain reliable results effectively. As given in Table 3, three types of load conditions: pure shear, y-directional uniaxial tension, and biaxial tension, are analyzed. The strain components of  x ,  y ,  xy 

T

for

various angles of anisotropy are calculated from the peridynamic simulation with Eq. (45) and compared to the respective theoretical values from Eq. (43). The strain components of orthotropic lamina under pure shear, y-directional tensile, and biaxial tensile loads are shown in Fig. 9. As shown in Fig. 9, the strains for those three kinds of load conditions are well predicted using the proposed peridynamic model with the maximum difference of 7.9 %. Since any types of static load condition can be expressed by the basic stress components  x ,  y , xy 

T

with the angle of material principal φ, the

proposed numerical model is thus capable of capturing the strain response of orthotropic lamina for any load conditions.

26

Table 3. Three types of load cases Cases

Stress states

Pure shear

 xy  1MPa,  x   y  0

y-directional uniaxial tension

 y  1MPa,  x   xy  0

Biaxial tension

 x   y  1MPa,  xy  0

24

12 x

8

 x

y

xy

PD

PD

PD

0

-8

PD

4

PD

y

Strain (× 10-3)

Strain (× 10-3)

16

8

xy PD

0 -4

0

15

30

45 

60

75

-8

90

0

15

30

(a)

45 

60

75

90

(b) 12

Strain (× 10-3)

6

0

 x

y

xy

PD

PD

PD

-6

-12

0

15

30

45 

60

75

90

(c) Fig. 9. Strain components of orthotropic plate under (a) pure shear, (b) y-directional tensile, and (c) biaxial tensile load.

28

4.2 Compact tension (CT) test The compact tension (CT) test is usually used to determine the mode Ι fracture toughness (critical strain energy release rate) of linear elastic isotropic metallic materials [25]. It has also been used to investigate the fracture characterizations of anisotropic materials [26,27]. In this section, the fracture behavior of an orthotropic material using the CT test is quantitatively analyzed using the present peridynamic model of orthotropic materials. uy(t)

2

1

y x 0.55 w

1.2 w

o

uy(t)

a w



1.25 w

Fig. 10. The compact tension (CT) specimen.

The sizes and load condition of CT specimen are shown in Fig. 10. According to [27], the relationship of the critical stress intensity factor, KIc, and critical applied load, Pc, is defined as K Ic 

Pc a f  1/2 Bw  w

(46)

where: 1/2

a a f    29.6    w  w

a  185.5    w

3/2

a  655.7    w

5/2

a  1017    w

7/2

a  638.9    w

9/2

(47)

a and w are defined in Fig. 10, and B is the thickness of the specimen. For orthotropic 29

materials, the relationships of KIc1 and KIc2 (crack on the plane normal to the 1- and 2-principal axes, respectively) with the critical energy released rate GIc1 and GIc2 can be expressed as [28]: 1/2

GIc1  K

2 Ic1

S11S22  S11 2S12  S66     , 2  S22 2S22 

(48)

1/2

GIc 2  K Ic2 2

S11S22  S22 2S12  S66     2  S11 2S11 

where Sij are the components of compliance matrix of orthotropic materials. The anisotropic material, cortical bone [27,29], is considered, and the material properties are reported in Table 4. The plane stress condition is investigated, and the constant thickness of B = 1 mm is used. In addition, the explicit time integration is utilized, and the uniform time step of 100 ns (nanosecond) is used in the whole example. The specimen is loaded with symmetrical and linearly-increasing displacement at a constant speed of 10 mm/s, as shown in Fig. 10. In this example, the fracture behavior of the CT test is first investigated with the proposed peridynamic model as a/w = 0.5 (Fig. 10). Then, the crack predictions of CT specimens with different pre-crack lengths or anisotropy (fiber) orientations are performed.

Table 4: Material properties of cortical bone [27,29] E1

E2

G12

(GPa)

(GPa)

(GPa)

22

11.3

5.4

v12

0.396

30

ρ

GIc1

GIc2

(kg/m3)

J/m2

J/m2

2000

1374

644

4.2.1 Simulations of CT test with a/w = 0.5 Similar to the CT test of cortical bone in [27], the typical load orientation condition (normal to the 1-principal material axis) is first considered as shown in Fig. 10, which leads to φ = 0°. The initial pre-crack length a0 = 20 mm and w = 40 mm are used, leading to a/w = 0.5. In the peridynamic model, the system is discretized into uniform grids, and the values of δ = 1 mm and 2 mm, and m = 4, 6 and 8 are performed for the sake of convergence study. The distributions of displacement in y-direction, damage (crack path), and released energy density at the typical state of CT crack propagation simulation are shown in Fig. 11, where the value of δ = 1 mm and m = 6 are used. As the y-directional displacement increases symmetrically (see Figs. 11(a) and (b)), the crack starts to grow from the location of the pre-crack tip (see Fig. 11(c)) and propagates along the pre-crack direction, paralleled to the 1-principal material axis (see Fig. 11(d)). Meanwhile, the released energy density appears at the pre-crack tip and follows with the crack path (see Figs. 11(e) and (f)). Typically, the contours of the strain energy density and released energy density around the location of the pre-crack tip before and after the crack growth are given in Fig. 12. As shown in Fig. 12, the strain energy density is first concentrated at pre-crack tip (see Fig.12(a)), and the value of the strain energy density increases with increasing displacement before crack starts to grow (see and compare Figs. 12(b) and (c)). Then, the crack begins to grow when the strain energy arrives at the critical value and the crack tip-concentrated location moves with the crack path (see Figs. 12(c) and (d)). Meantime, the released energy density appears before crack starts to grow (see Figs. 12(c) and (f)) and follows the crack tip-concentrated location of the strain energy density as crack propagates (see Figs. 12(d) and (g)).

31

(a) At 9 × 10-3 s

(b) At 20 × 10-3 s

(c) At 9 × 10-3 s

(d) At 20 × 10-3 s

(e) At 9 × 10-3 s

(f) At 20 × 10-3 s

Fig. 11. Distributions of y-direction displacement [(a) and (b)], damage/crack propagation [(c) and (d)], and released energy density [(e) and (f)] at the typical state of CT simulation.

32

(a) 6 × 10-3 s

(b) 7 × 10-3 s

(c) 8 × 10-3 s

(d) 9 × 10-3 s

(e) 7 × 10-3 s

(f) 8 × 10-3 s

(g) 9 × 10-3 s

Fig. 12. The contours of strain energy density [(a), (b), (c), and (d)] and released energy density [(e), (f), and (g)] around the location of crack tip of CT specimen.

The energy components changing with the increasing displacement loading are presented in Fig. 13. As shown in Fig. 13, the summation of Ue (the elastic strain energy), Uk (the kinetic energy), and WS (the incremental surface energy) is equal to We (the work done by external forces) during the CT simulation, which perfectly matches the energy balance 33

condition [30]:

dWe  dU e  dU k  dWS

(49)

12 Ue Uk Energy (× 10-3 J)

WS

8

Ue+Uk+WS We

4

0 0.00

0.05

0.10

0.15

0.20

uy (mm)

Fig. 13. Energy components during CT simulation with δ = 1 mm and m = 6.

For the convergence study results, the distributions of strain energy density and released energy density at the location of crack tip with a fixed horizon size δ = 1 mm and varying values of m = 4, 6 and 8 are shown in Fig. 14; while the distributions with a fixed value of m = 6 and varying values of δ = 2 mm and 1 mm are presented in Fig. 15. Overall, the strain energy and damage profiles for different values of m and δ have the similar pattern. As for the much denser grid size, i.e., larger values of m or smaller values of δ, the concentration effect is more obvious with smaller concentrated sizes and larger maximum values (Figs. 14 and 15). The load-displacement profiles are shown in Fig. 16. As shown in Fig. 16, the applied load P increases linearly with the increasing displacement load uy and drops with the oscillations as it arrives at the critical applied load. The plots of the load-displacement profiles for different values of m are nearly coincided and have the similar values of slope and critical applied load, which are reported in Table 5. Furthermore, the oscillations occur 34

when the cracks start to grow because of the explicit time integration considered in the transient analysis, and the oscillations of curves are weakened with the increasing value of m.

(a) m = 4

(b) m = 6

(c) m = 8

(e) m = 4

(f) m = 6

(g) m = 8

Fig. 14. Strain energy density [(a), (b) and (c)] and released energy density plots [(e), (f) and (g)] at 8 × 10-3 s with δ = 1 mm.

(a) δ = 2 mm

(b) δ = 1 mm

35

(c) δ = 2 mm

(d) δ = 1 mm

Fig. 15. Strain energy density [(a) and (b)] and released energy density plots [(c) and (d)] at 8 × 10-3 s with m = 6.

66 m=4 m=6 m=8

P (N)

44

22

0 0.00

0.05

0.10

0.15

0.20

u y (mm)

Fig. 16. Applied load vs displacement plots for the CT tests with δ = 1 mm.

The contours of the total incremental surface energy WS changing with the crack length a are shown in Fig. 17. As a whole, WS increases linearly with the crack length in the step-jump form. The calculated energy release rate is defined as:

GQ 

WS Ba

(50)

where B is the thickness of the 2D model. GQ is defined and compared to the input critical energy release rate in Table 4 to trace the released energy during crack propagation. The plots of the load-displacement profiles have the nearly same slope and converge to a straight line as

36

horizon decreases, and the values of slopes are also given in Table 5. The CT test simulation results for different values of δ and m are reported in Table 5. In Table 5, the CT compliance ue/P, the critical applied load Pc, and the calculated energy release rate GQ are given according to the peridynamic simulations, where ue is the displacement at the edge point o (see Fig. 10). The values of ue/P from the present peridynamic model closely match those from the numerical finite element method (FEM) model, with the maximum difference of 5.8%, where the CT compliance ue/P is an elastic parameter before crack starts to grow and it can be easily calculated using the symmetrical FEM model (see Fig. 10). The values of Pc are close to the standard values in Eq. (46) with the maximum difference of 3.9%. In addition, the maximum difference of GQ with the input critical energy release rate GIc2 in Table 4 is less than 6.0%. 8  = 2 mm  = 1 mm

WS (× 10-3 J)

6

4

2

0 30

33

36

39

42

a (mm)

Fig. 17. Incremental surface energy vs. crack length with m = 6. Table 5. Peridynamic simulation results for CT test with a/w = 0.5

?? (mm)

2

1 Classical

m ue / P (mm/kN)

4 1.87

6 1.85

8

4

1.83

1.82

37

6 1.80

8 1.79

1.78

Pc (N)

60.4

58.0

57.2

62.1

58.8

59.0

59.8

GQ (J/m2)

610

633

629

605

633

628

644

In summary, the proposed peridynamic model can accurately capture the fracture characteristics during the whole CT test (i.e., the compliance during loading before crack grows, the critical applied load when crack starts to grow, and the released energy balance during crack propagation).

4.2.2 CT test predictions with different pre-crack a/w In this section, a fixed value of w = 40 mm and varying values of a0 from 12 mm to 28 mm are considered, leading to different values of a/w. Also, the value of φ = 0° is considered, and the values of m = 6 and δ = 2 mm are utilized according to the convergence study given above. The compliance and critical applied load of CT specimens for different values of a/w are shown in Figs. 18 and 19, respectively. As shown in Fig. 18, the compliance of CT specimens increases exponentially with a/w, and the values of compliance from the present peridynamic model are close to those from the FEM model within a maximum difference of 4.7%. Meanwhile, as shown in Fig. 19, the critical applied load decreases with the increasing value of a/w, and the simulated values of critical applied load are close to those from Eq. (46) within a maximum difference of 3.9%.

38

8

u e / P (mm/KN)

6

PD Standard

4

2

0 0.3

0.4

0.5

0.6

0.7

a/w

Fig. 18. The compliance of CT specimens with different pre-crack a/w. 120 PD Standard

Pc (N)

80

40

0 0.3

0.4

0.5

0.6

0.7

a/w

Fig. 19. The critical applied load of CT specimens with different pre-crack a/w.

In addition, the quantitative values of parameters from the CT test simulations for different values of a/w are reported in Table 6. From Table 6, the values of calculated energy release rate GQ for different a/w closely match the input value of GIc2 (= 644 J/m2, see Table 4), with a maximum difference of 2.8%.

Table 6. Peridynamic simulation results for CT test with different a/w a/w

0.3

0.4

0.5

0.6

0.7

ue / P (mm/kN)

0.698

1.0

1.85

3.18

6.22

39

Pc (N)

99.2

78.5

58.0

40.7

25.8

GQ (J/m2)

632

630

633

631

626

4.2.3 CT test predictions with different material orientation φ Similar to the experiments in [27], different fiber orientations (anisotropy) φ = 15°, 30°, 45°, 60°, 75° and 90° are investigated for the CT specimen test (see Fig. 10). The initial pre-crack length a0 = 20 mm and w = 40 mm are used, and the values of m = 6 and δ = 2 mm are utilized. The crack path of the CT specimens with the different anisotropy (fiber) orientations are presented in Fig. 20. Similar to the results in [27], the crack grows approximately along the 1-principal material (or fiber) axis instead of the pre-crack direction, and the bifurcated fracture is obtained as φ = 90°. Since the released energy of a point varying with crack path orientation, the minimum absorbed energy path is naturally chosen as the crack propagation path, which is the 1-principal (fiber) direction in this example (see Table 4).

15°

(a) φ = 15°

30°

(b) φ = 30°

40

45°

(c) φ = 45°

60°

(e) φ = 60°

90°

75°

(f) φ = 75°

(g) φ = 90°

Fig. 20. Crack paths of the CT specimens with different anisotropy orientations φ.

4.3 Single edge-notched tension (SENT) test In this section, a single edge-notched tension (SENT) specimen test is further analyzed to validate the proposed orthotropic peridynamic model. The effect of critical stretches s2 and s1 on fracture analysis is investigated. The geometries and load condition of the SENT test are shown in Fig. 21, and the fixed pre-crack size a = 16 mm is utilized. The plane stress condition is considered, and the constant thickness of B = 1 mm is used. The orthotropic material of glass-epoxy unidirectional composite is considered, and the material constants are given in Table 7. The material fracture constant is not given here, and different values of critical stretches are discussed to investigate their effect on fracture analysis. For peridynamic simulation model of SENT specimen, the mesh sizes of m = 6 and δ = 2 mm are utilized.

41

u

2

y

1

x

80 mm

o

a

40 mm

u

Fig. 21. The single edge-notched tension (SENT) specimen. Table 7. Elastic properties of glass-epoxy composite [31] E1

E2

G12

(GPa)

(GPa)

(GPa)

40

8.3

4.1

v12

ρ (kg/m3)

0.26

1200

4.3.1 Elastic behavior of SENT with fiber orientation of φ = 0° The elastic behavior of SENT specimen in Fig. 21 is first analyzed, in which the anisotropy (fiber) orientation φ = 0° is considered. The fixed displacement load is uy = 4 × 10-3 m, and the ADR method is performed for the quasi-static analysis. The displacement and strain energy density by the peridynamic and FEM models are shown in Figs. 22 to 23, respectively. As shown in Figs. 22 to 23, the distributions of displacement and strain energy density by the peridynamic model are consistent with those predicted by FEM.

42

(a)

(b)

(c)

(d)

Fig. 22. Comparisons of displacements distributions in x-direction (a) FEM and (b) Peridynamics, and y-direction (c) FEM and (d) peridynamics of the SENT specimens with crack length a = 16 mm.

Fig. 23. Comparison of strain energy density distributions (a) FEM and (b) Peridynamics of the SENT specimens with crack length a = 16 mm.

43

4.3.2 Effect of s2 The crack growth prediction of SENT test is then investigated, and different given values of s1 and s2 are considered to evaluate their effect on fracture analysis. The uniform time step of 1.0 × 10-7 s for the explicit time integration is used. As shown in Fig. 21, the specimen is loaded with symmetrical and linearly-increasing displacement with a constant speed of 1.0 × 10-2 m/s. First, the fixed critical stretch values of s0 = 4.0 × 10-3 and s1 = 1.0 × 10-2 are utilized, and the varying values of s2 = 2.0 × 10-3, 4.0 × 10-3 and 6.0 × 10-3 are evaluated. The typical anisotropy (fiber) orientation φ = 0° is considered. The weights of broken bonds of d0, d1 and d2 before and after failure are shown in Fig. 24, where di is the point variable parameter that defined in Eq. (41). As shown in Fig. 24, with the case of the anisotropy (fiber) orientation φ = 0°, the crack grows along the crack original (1-principal axis) direction, where the d0 and d2 exist with the pre-crack (see Fig. 24(a) and (c)) and they develop as crack propagation (see Fig. 24(d) and (f)), indicating that the equivalent isotropic bonds  0 and the 2-axis additional bonds  2 are broken to form the crack. Meanwhile, d1 keeps unchanged (i.e., d1 = 0) during crack propagation (see Fig. 24(b) and (e)), demonstrating that the 1-axis additional bonds 1 are not broken.

(a) d0

(b) d1

44

(c) d2

(d) d0

(e) d1

(f) d2

Fig. 24. Broken weights of d0 [(a) and (d)], d1 [(b) and (e)] and d2 [(c) and (f)] of the SENT specimen with fiber orientation of φ = 0° before and after failure, respectively.

For the different values of s2 = 2.0 × 10-3, 4.0 × 10-3 and 6.0 × 10-3, the crack path are similar. However, the calculated energy released rates GQ from Eq. (50) are different as shown in Table 8, where the value of GQ decreases with the increasing value of s2. s2 has the negative relationship with the released energy, which may lead to energy inverse during crack propagation. Moreover, as aforementioned, the reasonable value of s2 = s0 is better to avoid such energy inverse (see the discussion after Eq. (39)). Table 8. Calculated energy released rates for SENT tests with different s2 s2

2 × 10-3

4 × 10-3

6 × 10-3

GQ (J/m2)

108.5

103.4

91.8

4.3.3 Effect of s1 Then, the fixed values of s0 = s2 = 4.0 × 10-3 are utilized, and the effect of s1 on crack growth prediction is evaluated. The different values of s1 = 4.0 × 10-3 and 1.0 × 10-2 are discussed, and the typical anisotropy (fiber) orientation φ = 45° is considered. The final weights of broken bonds of d0, d1 and d2 after failure are illustrated in Fig. 25. 45

As shown in Fig. 25, with the anisotropy (fiber) orientation φ = 45°, the pre-crack turns to grow along φ = 45° direction (i.e., still the 1-principal axis). Unlike the case of φ = 0°, d1 exists with the pre-crack and slightly appears along the crack path, but it is much smaller compared to d2 (see Fig. 25(b) and (e)). Specially, even for the same value of critical stretches (s1 = s0 = s2 = 4.0× 10-3), the crack still grows along the 1-axis direction (see Fig. 25(a) and (c)). Thus, not only the critical stretches but also the orthotropic constitutive properties impart effect on the crack path.

(a) d0 (s1 = 4.0 × 10-3)

(d) d0 (s1 = 1.0 × 10-2)

(b) d1 (s1 = 4.0 × 10-3)

(e) d1 (s1 = 1.0 × 10-2)

(c) d2 (s1 = 4.0 × 10-3)

(f) d2 (s1 = 1.0 × 10-2)

Fig. 25. Broken weights of d0 [(a) and (d)], d1 [(b) and (e)] and d2 [(c) and (f)] of the SENT specimens with φ = 45° for s1 = 4.0 × 10-3 and s1 = 1.0 × 10-2, respectively.

46

In addition, for different values of s1 = 4.0 × 10-3 and 1.0 × 10-2, the calculated energy released rates GQ have similar values as 113.9 and 113.3 J/m2. Thus, the value of s1 almost has no effect on released energy during crack propagation. In summary, the crack tends to grow along the 1-principle axis (fiber) direction even for the varying values of s1 and s2. The value of s2 has the negative relationship with the released energy as crack propagation, while the value of s1 has less effect on released energy. 5. Conclusions In this paper, a new state-based peridynamic model is proposed for quantitative fracture analysis of the two-dimensional (2D) orthotropic materials. The relationships of four peridynamic constants with the same number of orthotropic material parameters are for the first time established, and three types of critical stretches are given for fracture analysis of orthotropic materials. An orthotropic lamina under stress load and both the compact tension (CT) and single edge-notched tension (SENT) tests are quantitatively analyzed with the proposed model. The strain components of orthotropic lamina from peridynamic solutions converge to their respective theoretical results as the value of m increases or δ decreases. With the densest grid size, the maximum difference of numerical strain components with theoretical values is 6.1%, and the peridynamic model can still capture the strain response of orthotropic lamina for various load conditions. In the CT test case, the crack starts to grow from the location of the pre-crack tip and propagates parallel to the 1-principal axis. Also, the simulation model greatly captures the fracture characteristics (i.e., the crack path, the energy concentration, the displacement-load

47

relationship, and energy balance during crack propagation) during the whole CT test. For the different pre-crack lengths, the compliance and critical applied load of CT specimens are close to those from the numerical finite element models within the maximum difference of 4.7% and 3.9%, respectively. Meanwhile, with the different anisotropy (fiber) orientations, the crack of CT specimens chooses to grow along the 1-material principal (fiber) direction. In the SENT test, the crack also turns to grow along the 1-principal axis direction for both the cases of orientation φ = 0° and φ = 45°. The value of s2 has the negative relationship with the released energy density, while value of s1 almost has no effect on released energy. In summary, the proposed state-based peridynamic model for 2D orthotropic materials is capable of quantitatively capturing the elastic behaviors of orthotropic lamina and predicting fracture characteristics of CT and SENT specimens. Further, there are no extra failure criteria needed in the fracture analysis, and the proposed model can be effectively used to quantitatively analyze the fracture behavior of composite materials of anisotropic nature.

Acknowledgements The authors would like to acknowledge the partial financial support from the National Natural Science Foundation of China (NSFC Grant Nos. 51478265 and 51679136) to this study.

References [1] Kim JH, Paulino GH. Mixed-mode fracture of orthotropic functionally graded materials using finite elements and the modified crack closure method. Eng Fract Mech 2002;69:1557–86.

48

[2] Motamedi D, Mohammadi S. Dynamic crack propagation analysis of orthotropic media by the extended finite element method. Int J Fract 2010;161:21–39. [3] Aliabadi MH, Sollero P. Crack growth analysis in homogeneous orthotropic laminates. Compos Sci Technol 1998;58:1697–703. [4] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 2000;48:175–209. [5] Silling SA, Epton M, Weckner O, Xu J, Askari E. Peridynamic states and constitutive modeling. J Elast 2007;88:151-184. [6] Zhang H, Qiao P. An extended state-based peridynamic model for damage growth prediction of bimaterial structures under thermomechanical loading. Eng Fract Mech 2017;189:81–97. [7] Ha YD, Bobaru F. Characteristics of dynamic brittle fracture captured with peridynamics. Eng Fract Mech 2011;78:1156–68. [8] Huang D, Lu G, Qiao P. An improved peridynamic approach for quasi-static elastic deformation and brittle fracture analysis. Eng Fract Mech 2015;94–95:111–22. [9] Askari E, Xu J, Silling S. Peridynamic analysis of damage and failure in composites. In: 44th AIAA/ASME/ASCE/AHS/ASC aerospace sciences meeting and exhibit, 9–12 January 2006, Reno, NV. AIAA 2006-88. [10] Gerstle W, Sau N, Silling S. Peridynamic modeling of plain and reinforced concrete structures. In: SMiRT18. Int conf structural mech reactor technol; 2005. [11] Xu J, Askari A, Weckner O, Razi H, Silling S. Damage and failure analysis of composite laminates under biaxial loads. In: 48th AIAA/ASME/ASCE/AHS/ASC structures,

49

structural dynamics, and materials conference, 23–26 April 2007, Honolulu, HI. AIAA 2007-2315. [12] Oterkus E, Madenci E, Weckner O, Silling S, Bogert P, Tessler A. Combined finite element and peridynamic analyses for predicting failure in a stiffened composite curved panel with a central slot. Compos Struct 2012;94:839–50. [13] Hu W, Ha YD, Bobaru F. Peridynamic model for dynamic fracture in unidirectional fiber-reinforced composites. Comput Methods Appl Mech Eng 2012;217–220:247–61. [14] Ghajari M, Iannucci L, Curtis P. A peridynamic material model for the analysis of dynamic crack propagation in orthotropic media. Comput Methods Appl Mech Eng 2014;276:431–52. [15] Madenci E, Oterkus E. Peridynamic Theory and Its Applications. Springer, New York; 2014. [16] Silling SA. Linearized theory of peridynamic states. J Elast 2010;99:85–111. [17] Zhang H, Qiao P. A state-based peridynamic model for quantitative fracture analysis. Int J Fract 2018;211:217–35. [18] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Comput Struct 2005;83:1526–35. [19] Bobaru F, Yang M, Alves LF, Silling SA, Askari E, Xu J. Convergence, adaptive refinement, and scaling in 1D peridynamics. Int J Numer Methods Eng 2009;77:852–877 [20] Le Q V, Bobaru F. Surface corrections for peridynamic models in elasticity and fracture. Comput Mech 2017:1–20. [21] Bobaru F, Ha YD. Adaptive Refinement and Multiscale Modeling in 2D Peridynamics.

50

Int J Multiscale Comput Eng 2011;9:635–60. [22] Hu YL, Madenci E. Bond-based peridynamic modeling of composite laminates with arbitrary fiber orientation and stacking sequence. Compos Struct 2016;153:139–75. [23] Kilic B, Madenci E. An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory. Theor Appl Fract Mech 2010;53:194–204. [24] Sarego G, Le QV, Bobaru F, Zaccariotto M, Galvanetto U. Linearized state-based peridynamics for 2-D problems. Int J Numer Methods Eng 2016;108:1174–97. [25] ASTM E399-12e3. Standard test method for linear elastic plane strain fracture toughness KIC of metallic material; 2013. [26] Norman TL, Vashishth D, Burr DB. Fracture toughness of human bone under tension. J Biomech 1995;28. [27] Behiri JC, Bonfield W. Orientation Dependence of the Fracture Mechanics of Cortical Bone. J Biomech 1989;22. [28] Sih GC, Paris PC, Irwin GR. On cracks in rectilinearly anisotropic bodies. Int J Fract Mech 1965;1:189–203. [29] Lang SB. Ultrasonic Method for Measuring Elastic Coefficients of Bone and Results on Fresh and Dried Bovine Bones. IEEE Trans Biomed Eng 1970;17:101–5. [30] Sun CT. Fracture Mechanics. Academic Press,Waltham, MA; 2012. [31] Parhizgar S, Zachary LW, Sun CT. Application of the principles of linear fracture mechanics to the composite materials. Int J Fract 1982;20:3-15.

51

A state-based peridynamic model for quantitative elastic and fracture analysis of orthotropic materials Heng Zhanga, Pizhong Qiaoa,b, a

State Key Laboratory of Ocean Engineering, Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, PR China b

Department of Civil and Environmental Engineering, Washington State University, Sloan Hall 117, Pullman, WA 99164-2910, USA

Highlights: -

A new state-based peridynamic model is proposed for quantitative elastic and fracture analyses of the orthotropic materials.

-

The relationships between four independent material parameters of two-dimensional (2D) orthotropic materials and the same number of peridynamic constants are for the first time established.

-

Three types of critical stretches are uniquely defined for quantitative fracture analysis of orthotropic materials.

-

The model is capable of quantitatively capturing fracture characteristics of orthotropic materials.



Corresponding author. E-mail addresses: [email protected] and [email protected] (P. Qiao). 52