Peridynamics for bending of beams and plates with transverse shear deformation

International Journal of Solids and Structures xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr

Peridynamics for bending of beams and plates with transverse shear deformation C. Diyaroglu a, E. Oterkus a, S. Oterkus b, E. Madenci b,⇑ a b

Department of Naval Architecture, Ocean and Marine Engineering, University of Strathclyde, Glasgow, United Kingdom Department of Aerospace and Mechanical Engineering, University of Arizona, Tucson, United States

a r t i c l e

i n f o

Article history: Received 2 October 2014 Received in revised form 12 April 2015 Available online xxxx Keywords: Peridynamics Timoshenko beam Mindlin plate Transverse shear deformation Dispersion relationships

a b s t r a c t Progressive failure analysis of structures is still a major challenge. There exist various predictive techniques to tackle this challenge by using both classical (local) and nonlocal theories. Peridynamic (PD) theory (nonlocal) is very suitable for this challenge, but computationally costly with respect to the finite element method. When analyzing complex structures, it is necessary to utilize structural idealizations to make the computations feasible. Therefore, this study presents the PD equations of motions for structural idealizations as beams and plates while accounting for transverse shear deformation. Also, their PD dispersion relations are presented and compared with those of classical theory. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Peridynamic (PD) theory was originally introduced for the solution of deformation field equations (Silling, 2000) without any structural idealizations. It satisfies all the fundamental balance laws of classical (local) continuum mechanics; however, it is different in the sense that it is a nonlocal continuum theory and it introduces an internal length parameter into the field equations. This internal length parameter defines the association among the material points within a finite distance through micropotentials. Removal of micropotentials between the material points allows damage initiation and growth through a single critical failure parameter regardless of the mixed-mode loading conditions. The creation of a new (crack) surface is based on a local damage measure. The local damage is defined as the ratio of broken interactions to the total number of interactions at a material point. Finite element analysis (FEA) with traditional elements suffers from the following shortcomings: (1) the interface between dissimilar materials is assumed to have zero thickness without any specific material properties; however, it presents a weak link and it is usually the location of failure. Therefore, it fails to appropriately model the interface between dissimilar materials. (2) Failure is a dynamic process, and it requires remeshing. It is computationally costly, and the crack growth is guided based on the linear elastic fracture mechanics (LEFM) concepts. It breaks down when multiple complex crack growth patterns develop. (3) Stress ⇑ Corresponding author.

and strain fields are discontinuous, and mesh refinement does not necessarily ensure accurate stress fields near geometric and material discontinuities. (4) Finally, crack nucleation is not resolved. The analysis always requires a pre-existing crack. In order to remedy or remove these shortcomings, Cohesive Zone Elements (CZE) and eXtended Finite Elements (XFEM) were developed; however, CZE requires a priori knowledge of the crack path. In a complex analysis, it is not practical and the results are dependent on the mesh (structured or unstructured). Furthermore, the results are sensitive to the strength parameters in the traction–separation law of the cohesive zone model. Determination of these parameters poses additional uncertainties. Although XFEM removed such uncertainties, it still requires an external criteria for crack propagation. Thus, the results depend on the criteria employed in the analysis. It also breaks down when multiple complex crack growth patterns develop. The PD theory overcomes the weaknesses of the existing methods, and it is capable of identifying all of the failure modes without simplifying assumptions. The PD methodology effectively predicts complex failure in complex structures under general loading conditions. Damage is inherently calculated in a PD analysis without special procedures, making progressive failure analysis more practical. An extensive literature survey on PD is given in a recently published textbook by Madenci and Oterkus (2014). A comparison study between peridynamics, CZE, and XFEM techniques is given by Agwai et al. (2011). They showed that the crack speeds obtained from all three approaches are on the same order; however, the

http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040 0020-7683/Ó 2015 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

2

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

fracture paths obtained by using peridynamics are closer to experimental results with respect to other two techniques. Another advantage of PD is its length-scale parameter, which does not exist in classical continuum mechanics. Such a length-scale parameter gives PD a nonlocal character. Hence, it allows the capture of physical phenomena not only at the macro-scale, but also at various other scales. This characteristic can be established through the PD dispersion relations. The classical theory is only valid for a special case of a long wavelength limit; however, the PD shows dispersion behavior similar to that observed in real materials. Hence, it is proven to be acceptable to perform multi-scale analysis simulations. Although peridynamics is a powerful technique in failure analysis and has an internal length scale, it is usually computationally more expensive, especially with respect to finite element analysis. The computational time can be significantly reduced by using parallel computing either by using a CPU (Central Processing Unit) and/or GPU – based (Graphics Processing Unit) architecture for which PD equations of motion are very suitable. However, modeling very large and detailed structures such as aerospace and marine vehicles can still be computationally demanding. Hence, in such cases it is necessary to reduce computational time through structural idealization. Taylor and Steigmann (2013) proposed a peridynamic plate model based on bond-based formulation by using an asymptotic analysis. Their formulation is capable of capturing out-of-plane deformations for thin plates. Moreover, O’Grady and Foster (2014a,b) developed a non-ordinary state-based peridynamic model for Euler–Bernoulli beam and Kirchhoff–Love plate formulations by disregarding the transverse shear deformations. Therefore, the focus of this study is to present a new PD formulation for thin or thick beams and plates while accounting for transverse shear deformation based on an original (bond-based) PD formulation. Moreover, PD dispersion relations are obtained and compared against those from classical theory. The following sections present the PD kinematics for a Timoshenko beam and a Mindlin plate, and the corresponding PD equations of motion as well as the PD material parameters. They also describe the procedure to determine the surface correction factors for these parameters and the application of the boundary conditions and determination of the critical curvature and critical shear angle in terms of the fracture mechanics parameters. Finally, the corresponding dispersion relations are derived and compared with the classical theory. The numerical results establish the validity of the present formulation by considering simple benchmark problems.

2. Peridynamic kinematics At any instant of time, every point in the beam or plate denotes the out-of-plane deflection and rotations of a material particle, and these infinitely many material points (particles) constitute the beam or the plate. In the undeformed state of the body, each material point is identified by its coordinates, xðkÞ with ðk ¼ 1; 2; . . . ; 1Þ, and is associated with an incremental volume, V ðkÞ , and a mass density of qðxðkÞ Þ: According to the PD theory introduced by Silling (2000), the motion of a body is analyzed by considering the pair-wise interaction between material points xðkÞ and xðjÞ . The interaction between the material points is prescribed through a micropotential that depends on the deformation and constitutive properties of the material. Also, a material point is only influenced by the other material points within a neighborhood defined by its horizon, d. The micropotentials are zero for material points outside its horizon. Each material point can be subjected to prescribed body loads, displacement, or velocity, resulting in motion and deformation.

2.1. Beam kinematics As shown in Fig. 1, the transverse shear angles, uðjÞ and uðkÞ , of material points j and k can be expressed as

uðjÞ ¼ uðkÞ ¼

  wðjÞ  wðkÞ  /ðjÞ sgnðxðjÞ  xðkÞ Þ nðjÞðkÞ 

wðjÞ  wðkÞ  /ðkÞ sgnðxðjÞ  xðkÞ Þ nðjÞðkÞ

ð1aÞ

 ð1bÞ

in which wðjÞ , /ðjÞ and wðkÞ , /ðkÞ represent the out-of-plane deflection and rotation of material points j and k, respectively. The distance   between the material points j and k is specified as nðjÞðkÞ ¼ xðjÞ  xðkÞ . Considering the material point k as the point of interest, the transverse shear angle, uðkÞðjÞ , arising from the interaction between material points j and k can be defined as the average of the transverse shear angles at these material points in the form

uðkÞðjÞ ¼

  wðjÞ  wðkÞ /ðjÞ þ /ðkÞ  sgnðxðjÞ  xðkÞ Þ nðjÞðkÞ 2

ð2Þ

The curvature between the material points j and k can be defined as

jðkÞðjÞ ¼

  /ðjÞ  /ðkÞ nðjÞðkÞ

ð3Þ

When considering the material point j as the point of interest, the transverse shear angle and curvature for the interaction between the material points j and k become

uðjÞðkÞ ¼

    /ðkÞ þ /ðjÞ wðkÞ  wðjÞ sgnðxðjÞ  xðkÞ Þ or   nðjÞðkÞ 2

uðjÞðkÞ ¼ uðkÞðjÞ

ð4aÞ

and

jðjÞðkÞ ¼

  /ðkÞ  /ðjÞ or nðjÞðkÞ

jðjÞðkÞ ¼ jðkÞðjÞ

ð4bÞ

2.2. Plate kinematics As illustrated in Fig. 2, /ðjÞ and /ðkÞ represent the rotations with respect to the line of action between the material points j and k. Considering the material point k as the point of interest, the curvature, jðkÞðjÞ , with respect to the line of action between the material points j and k can be defined as

jðkÞðjÞ ¼

/ðjÞ  /ðkÞ nðjÞðkÞ

ð5Þ

Fig. 1. Original and deformed configurations of a Timoshenko beam.

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

3

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

3. Peridynamic equations of motion The PD equations of motion at material point k can be derived by applying the principle of virtual work

Z

t1

d

ðT  UÞdt ¼ 0

ð12Þ

t0

where T and U represent the total kinetic and potential energies in the beam or plate. This principle is satisfied by solving for the Lagrange equation

  d @L @L  ¼0 dt @ q_ ðkÞ @qðkÞ

Fig. 2. Original and deformed configurations of a Mindlin plate.

ð13Þ

Through coordinate transformation, the rotations and curvature with respect to the line of action between the material points j and k can be decomposed as

where the vector qðkÞ includes the independent field variables (out-of-plane deflection and rotations), and the Lagrangian L is defined as

/ðjÞ ¼ /xðjÞ cos h þ /yðjÞ sin h

ð6aÞ

L¼T U

/ðkÞ ¼ /xðkÞ cos h þ /yðkÞ sin h

ð6bÞ

3.1. Beam equations of motion

and

jðkÞðjÞ

!   /xðjÞ  /xðkÞ /yðjÞ  /yðkÞ 2 sin h cos2 h þ ¼ xðjÞ  xðkÞ yðjÞ  yðkÞ

ð7Þ

in which xðjÞ  xðkÞ ¼ nðjÞðkÞ cos h and yðjÞ  yðkÞ ¼ nðjÞðkÞ sin h, with nðjÞðkÞ representing the distance between material points j and k. The slope with respect to the line of action between the material points j and k can be expressed as

#ðkÞðjÞ ¼

wðjÞ  wðkÞ nðjÞðkÞ

ð8Þ

in which wðjÞ and wðkÞ represent the out-of-plane deflections at material points j and k. As sketched in Fig. 2, the transverse shear angles at material points j and k can be expressed as

uðjÞ ¼ #ðkÞðjÞ  /ðjÞ

ð9aÞ

uðkÞ ¼ #ðkÞðjÞ  /ðkÞ

ð9bÞ

Considering the material point k as the point of interest, the transverse shear angle, uðkÞðjÞ , between material points j and k can be defined as the average of the shear angles at these material points in the form

uðkÞðjÞ ¼

ð11aÞ

and

/ðkÞ  /ðjÞ nðjÞðkÞ

It is worth noting that uðjÞðkÞ ¼ uðkÞðjÞ and

ð15Þ

in which V ðkÞ represents the infinitesimally small incremental volume of material point k and the dot () above a parameter denotes differentiation with respect to time. The parameters q, I, and A correspond to mass density of the material, the moment of inertia, and the cross sectional area of the beam, respectively. The total potential energy of the system can be obtained by ~ ðkÞðjÞ ðjðkÞðjÞ Þ and w ^ ðkÞðjÞ ðuðkÞðjÞ Þ, summing the micropotentials, w between material points arising from bending and transverse shear deformation

( ) 1 1 X

1X 1 ~ ~ ðkÞðjÞ ðjðkÞðjÞ Þ þ w ~ ðjÞðkÞ ðjðjÞðkÞ Þ V ðjÞ  bðkÞ /ðkÞ V ðkÞ U¼ w 2 j¼1 2 k¼1 ( ) 1 1 i X 1X 1h ^ ^ ^ þ wðkÞðjÞ ðuðkÞðjÞ Þ þ wðjÞðkÞ ðuðjÞðkÞ Þ V ðjÞ  bðkÞ wðkÞ V ðkÞ 2 j¼1 2 k¼1

d @L @L  ¼0 _ ðkÞ @wðkÞ dt @ w

When considering the material point j as the point of interest, the transverse shear angle and curvature for the interaction between the material points j and k become

jðjÞðkÞ ¼

 1  1X I  q w_ 2ðkÞ þ /_ 2ðkÞ V ðkÞ 2 k¼1 A

~ and b ^ represent the body moment and body force at in which b ðkÞ ðkÞ material point k. The independent variables are out-of-plane deflection and rotation of the material point, wðkÞ and /ðkÞ . Hence, the resulting Euler–Lagrange equations can be expressed as

wðjÞ  wðkÞ nðjÞðkÞ     /xðjÞ cos h þ /yðjÞ sin h þ /xðkÞ cos h þ /yðkÞ sin h  2

  /ðkÞ þ /ðjÞ wðkÞ  wðjÞ   nðjÞðkÞ 2

T¼

ð10aÞ

ð10bÞ

uðjÞðkÞ ¼

The total kinetic energy of the system due to bending and transverse shear deformation can be written as

ð16Þ

wðjÞ  wðkÞ /ðjÞ þ /ðkÞ  nðjÞðkÞ 2

or

uðkÞðjÞ ¼

ð14Þ

and

d @L @L  ¼0 dt @ /_ ðkÞ @/ðkÞ

ð17bÞ

Using the Lagrangian definition L ¼ T  U and performing differentiation yield the following equations of motion

€ ðkÞ þ qw ð11bÞ

jðjÞðkÞ ¼ jðkÞðjÞ .

ð17aÞ

  1 X @ uðkÞðjÞ @ uðjÞðkÞ 1 ^ ¼0 V ðjÞ  b nðjÞðkÞ ^f ðkÞðjÞ þ nðjÞðkÞ ^f ðjÞðkÞ ðkÞ 2 @ðwðkÞ Þ @ðwðkÞ Þ j¼1 ð18aÞ

and

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

4

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

" ~f

# @ðjðkÞðjÞ Þ ~ @ðjðjÞðkÞ Þ V ðjÞ /ðkÞ þ n þ f ðjÞðkÞ ðkÞðjÞ A 2 ðjÞðkÞ @ð/ðkÞ Þ @ð/ðkÞ Þ j¼1 " # 1 X @ uðkÞðjÞ ^ @ uðjÞðkÞ 1 ~ ¼0 þ nðjÞðkÞ ^f ðkÞðjÞ þ f ðjÞðkÞ V ðjÞ  b ðkÞ 2 @ð/ðkÞ Þ @ð/ðkÞ Þ j¼1 1 X 1

qI €

3.2. Plate equations of motion

ð18bÞ

The total kinetic energy of the system due to bending and transverse shear deformation can be expressed as

T¼

in which ^f ðkÞðjÞ , ^f ðjÞðkÞ , ~f ðkÞðjÞ , and ~f ðjÞðkÞ are defined as

^f ðkÞðjÞ ¼

1

^ ðkÞðjÞ ðuðkÞðjÞ Þ @w

nðjÞðkÞ

@ uðkÞðjÞ

;

^f ðjÞðkÞ ¼

1

^ ðjÞðkÞ ðuðjÞðkÞ Þ @w

nðjÞðkÞ

@ uðjÞðkÞ

ð19a; bÞ

and

~f ðkÞðjÞ ¼

1 nðjÞðkÞ

~ ðkÞðjÞ ðjðkÞðjÞ Þ @w ; @ðjðkÞðjÞ Þ

~f ðjÞðkÞ ¼

1 nðjÞðkÞ

~ ðjÞðkÞ ðjðjÞðkÞ Þ @w @ðjðjÞðkÞ Þ

ð20a; bÞ

ð21a; bÞ

and

~f and ~f ðjÞðkÞ ¼ cb ðjðjÞðkÞ Þ ðkÞðjÞ ¼ c b ðjðkÞðjÞ Þ

ð22a; bÞ

in which cs and cb are the peridynamic material parameters associated with the transverse shear deformation and bending of the beam, respectively.Invoking Eqs. (2)–(4) and substituting for the peridynamic forces from Eqs. (21) and (22) in Eq. (18) result in

€ ðkÞ ¼ cs qw

Invoking the representation of the in-plane displacement components in terms of rotations as uðkÞ ¼ z/xðkÞ and v ðkÞ ¼ z/yðkÞ and performing integration through the thickness of the plate result in

! Z h=2 h 1 i 1 X 2 2 _2 2 _2 _ T¼ q w þ z /xðkÞ þ z /yðkÞ dz AðkÞ 2 k¼1 h=2 ðkÞ

 1  X wðjÞ  wðkÞ /ðjÞ þ /ðkÞ ^  sgnðxðjÞ  xðkÞ Þ V ðjÞ þ b ðkÞ nðjÞðkÞ 2 j¼1

! 2 2 1 1 X h _2 h _2 2 _ AðkÞ / þ / wðkÞ þ T ¼ hq 2 k¼1 12 xðkÞ 12 yðkÞ

and

A

/ðkÞ ¼ cb

1 X /ðjÞ  /ðkÞ

1  wðjÞ  wðkÞ 1 X V ðjÞ þ cs sgnðxðjÞ  xðkÞ Þ 2 j¼1 nðjÞðkÞ

nðjÞðkÞ  /ðjÞ þ /ðkÞ ~ nðjÞðkÞ V ðjÞ þ b  ðkÞ 2 j¼1

ð23bÞ

The peridynamic interactions are disregarded beyond the horizon. Therefore, these equations can be rewritten in an integral form over the domain of interactions as

Z 

wðx0 ; tÞ  wðx; tÞ n H  /ðx0 ; tÞ þ /ðx; tÞ 0 ^ tÞ sgnðx0  xÞ dV þ bðx;  2

€ tÞ ¼ cs qwðx;

ð24aÞ

ð28Þ

2k G Ad2

and cb ¼

2E I d2 A2

þ

1 kG 4 A

d @L @L  ¼0 _ ðkÞ @wðkÞ dt @ w

ð29aÞ

and

ð29bÞ

Using the Lagrangian definition L ¼ T  U, performing the differentiation while invoking Eqs. (6), (7), (10) and (11), and substituting from Eqs. (19)–(22) yield the following equations of motion

€ ðkÞ ¼ cs qhw

ð24bÞ

q

As derived in Appendix A, the PD material constants cs and cb can be expressed in terms of the shear and Young’s moduli, G and E, as

cs ¼

~ ^ in which b aðkÞ and bðkÞ represent the resultant body moment and body force at material point k. The independent variables are out-of-plane deflection and rotations of the material point, wðkÞ and /aðkÞ . Hence, the resulting Euler–Lagrange equations can be expressed as

1  X wðjÞ  wðkÞ /xðjÞ þ /xðkÞ  cos h nðjÞðkÞ 2 j¼1  /yðjÞ þ /yðkÞ ^  sin h V ðjÞ þ b ðkÞ 2

qI €

ð27bÞ

( ) 1 1 ~ X

b 1X 1 aðkÞ ~ ~ V ðkÞ U¼ / wðkÞðjÞ ðjðkÞðjÞ Þ þ wðjÞðkÞ ðjðjÞðkÞ Þ V ðjÞ  2 j¼1 2 h aðkÞ k¼1 ( ) 1 1 i ^ X b 1X 1h ^ ðkÞðjÞ ðuðkÞðjÞ Þ þ w ^ ðjÞðkÞ ðuðjÞðkÞ Þ V ðjÞ  ðkÞ wðkÞ V ðkÞ þ w 2 j¼1 2 h k¼1

d @L @L  ¼ 0 with a ¼ x; y dt @ /_ aðkÞ @/aðkÞ

and

 Z  /ðx0 ; tÞ  /ðx; tÞ 0 dV /ðx; tÞ ¼ cb A n H  Z  1 wðx0 ; tÞ  wðx; tÞ þ cs sgnðx0  xÞ 2 n H   /ðx0 ; tÞ þ /ðx; tÞ 0 ~ tÞ n dV þ bðx;  2

ð27aÞ

in which AðkÞ represents the infinitesimally small incremental area of each material point and h represents the thickness of the plate. The total potential energy of the plate can be obtained by sum~ ðkÞðjÞ ðjðkÞðjÞ Þ and w ^ ðkÞðjÞ ðuðkÞðjÞ Þ, between ming the micropotentials, w material points arising from bending and transverse shear deformation

ð23aÞ

qI €

ð26Þ

or

They represent the peridynamic interaction forces between material points j and k arising from transverse shear deformation and bending. For a linear material behavior, they can also be defined in the form

^f and ^f ðjÞðkÞ ¼ cs ðuðjÞðkÞ Þ ðkÞðjÞ ¼ c s ðuðkÞðjÞ Þ

1 h i 1 X _ 2ðkÞ V ðkÞ q u_ 2ðkÞ þ v_ 2ðkÞ þ w 2 k¼1

ð25a; bÞ

in which k is the shear correction factor, and is equal to 5=6 for rectangular cross sections.

ð30aÞ

    3 1  X /xðjÞ  /xðkÞ /yðjÞ  /yðkÞ h € /xðkÞ ¼ cb cos h þ sin h cos hV ðjÞ 12 nðjÞðkÞ nðjÞðkÞ j¼1  1 wðjÞ  wðkÞ /xðjÞ þ /xðkÞ 1 X þ cs nðjÞðkÞ  cos h 2 j¼1 nðjÞðkÞ 2  /yðjÞ þ /yðkÞ ~  ð30bÞ sin h cos hV ðjÞ þ b xðkÞ 2

and

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

5

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

q

    3 1  X /xðjÞ  /xðkÞ /yðjÞ  /yðkÞ h € /yðkÞ ¼ cb cos h þ sin h sin h V ðjÞ 12 nðjÞðkÞ nðjÞðkÞ j¼1  1 wðjÞ  wðkÞ /xðjÞ þ /xðkÞ 1 X þ cs nðjÞðkÞ  cos h 2 j¼1 nðjÞðkÞ 2  /yðjÞ þ /yðkÞ ~ ð30cÞ  sin h sin h V ðjÞ þ b yðkÞ 2

Note that the peridynamic interactions exist only within the horizon of material points. As derived in Appendix A, the PD material constants cs and cb can be expressed in terms of the shear and Young’s moduli, G and E, as

cs ¼

9E

2

k 4pd3

2

and cb ¼

E 3h 27 2 þ k pd 4d2 80

! ð31a; bÞ

this layer has a thickness defined by a single spacing between the material points. 5. Peridynamic dispersion relations Peridynamic dispersion relations are compared against those of the classical Timoshenko and Mindlin plate theories. In the derivation of these dispersion relations, the wave number, the wave frequency and phase velocity of the wave are denoted by j, x, and t, respectively. The relationship among these parameters is x ¼ j t. The compressional and shear wave speeds are defined by t2c ¼ E=q and t2s ¼ G=q, respectively, where G and E are the shear and Young’s moduli of the material. The wave number is related to the half-wavelength, C, by the relationship j ¼ p=C. 5.1. Beam dispersion relations

2

in which k is the shear correction factor that can be chosen based 2

on the frequency of the lowest thickness shear mode as k ¼ p2 =12, and m ¼ 1=3. The PD material parameters cb and cs are determined for a material point whose horizon is completely embedded in the material. For these material points, both classical and peridynamic Strain Energy Densities (SED) are equivalent. However, the PD material parameters cb and cs require correction if the material point is close to boundaries. The correction of the material parameters is achieved by numerically integrating the SED at each material point inside the body for simple loading conditions and comparing them to their counterparts obtained from classical solutions. SED is composed of bending and transverse shear deformation. Therefore, the parameter cb is corrected by gb based on the SED due to the bending, and cs by gs due to the transverse shear deformation. Their derivation is explicitly described in Appendix B.

Dispersion relations are determined by considering a wave propagating in the x-direction. Therefore, wave solutions for material points located at x and x0 can be expressed as 0

wðx; tÞ ¼ w0 eiðjxxtÞ

and wðx0 ; tÞ ¼ w0 eiðjxxtþsgnðx xÞjnÞ

/ðx; tÞ ¼ /0 eiðjxxtÞ

and /ðx0 ; tÞ ¼ /0 eiðjxxtþsgnðx xÞjnÞ

0

ð32aÞ ð32bÞ

in which jn is the phase difference between material points located at x and x0 , and /0 and w0 represent the amplitudes of these waves. Substituting these wave solutions into the PD equations of motion given in Eq. (24) leads to a homogeneous set of equations for /0 and w0 . For a nontrivial solution to exist, the determinant of their coefficient matrix must vanish, resulting in the wave dispersion relation as

  2Ac B  qx2 s 2    Acs B3

   ¼0 1 I 2 2Acb B2 þ 2 Acs B1  q A x Acs B3

ð33Þ

4. Boundary conditions Unlike the local theory, the boundary conditions are imposed through a nonzero volume of fictitious boundary layers. This necessity arises because the PD field equations do not contain any spatial derivatives; therefore, constraint conditions are, in general, not necessary for the solution of an integro-differential equation. However, such conditions can be imposed by prescribing constraints through a fictitious boundary layer. In order to apply a displacement or rotation constraint, a fictitious boundary layer, Rc , is introduced outside the actual material, as shown in Fig. 3. The size of this layer is equivalent to the horizon. An external load, such as a moment or a transverse load, can be applied in the form of body loads through a layer within the actual material, R‘ . As explained by Madenci and Oterkus (2014),

(a)

where the terms Bi with i ¼ 1; 2; 3 are explicitly given in Appendix C; they are dependent on the phase difference and the horizon. As explained in Appendix C, for long wavelength (or small wave number, j), the resulting wave dispersion relation is the same as that of the classical theory (Reis, 1978; Amirkulova, 2011). As expected, both theories yield the same relationship for long wavelength. For specified values of E ¼ 200 GPa, q ¼ 7850 kg=m3 , k ¼ 5=6, h ¼ 107 m, m ¼ 0:3, and finite horizon size, d ¼ 108 m, the evaluation of the determinant without any simplification leads to the variation of the wave frequency, x, as a function of the wave number, j, for the first and second modes as shown in Fig. 4. Although both modes involve a combination of transverse and angular displacements, the first mode is dominated by transverse motion and the second mode is dominated by angular motion (Reis,

(b)

Fig. 3. Application of boundary conditions in peridynamics: (a) beam and (b) plate.

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

6

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Fig. 4. Comparison of PD and CT wave frequency dispersions: (a) first mode and (b) second mode.

Fig. 5. Comparison of PD and CT wave speed dispersions: (a) first mode and (b) second mode.

1978). For both modes, the PD wave dispersions level off as the wave number increases which is a well-known behavior observed in experimental studies (Eringen, 1972; Weckner and Silling, 2011). However, the wave frequency always increases linearly according to the classical theory (CT). pffiffiffi The variation of the normalized phase speed, t=ðts kÞ as a function of wave number, j, for the first mode is shown in Fig. 5a. As shown in Fig. 5a, both PD and CT predict zero speed in the limit as j approaches zero. The phase speed of the CT reaches a constant value close to the shear wave speed of a bar for short wavelengths (or relatively large wave numbers). For the second mode, both theories predict comparable results for the long wavelengths (or relatively small wave numbers) while classical phase speed reaches the compressional wave speed, tc , of a bar for short wavelengths (or relatively large wave numbers), as depicted in Fig. 5b. However, the phase speed decreases as the wave number increases according to the PD theory as observed in real materials. 5.2. Plate dispersion relations As in the case of a beam, the dispersion relations for a plate can be obtained by considering a wave propagating in the x-direction. Therefore, wave solutions for material points located at x and x0 can be expressed as

wðx; tÞ ¼ w0 eiðjxxtÞ /x ðx; tÞ ¼ /x0 eiðjxxtÞ

and wðx0 ; tÞ ¼ w0 eiðjxxtþjn cosðhÞÞ and /x ðx0 ; tÞ ¼ /x0 eiðjxxtþjn cosðhÞÞ

ð34aÞ ð34bÞ

/y ðx; tÞ ¼ /y0 eiðjxxtÞ

and /y ðx0 ; tÞ ¼ /y0 eiðjxxtþjn cosðhÞÞ

ð34cÞ

in which jn cosðhÞ is the phase difference between material points located at x and x0 , and w0 , /x0 and /y0 represent the amplitudes of these waves. Substituting these wave solutions into the PD equations of motion given in Eq. (30) leads to a homogeneous set of equations for /x0 , /y0 , and w0 . For a nontrivial solution to exist, the determinant of their coefficient matrix must vanish, resulting in the wave dispersion relation as

    cs M 1 þ qhx2 cs M2 cs M 3     3 h 1 ¼0  cM cb M 4  12 cs M 7 þ q12 x2 cb M 5  12 cs M 8 2 s 6     3 qh 1 1 2  1 cs M 10 cb M 5  2 cs M 8 cb M9  2 cs M11 þ 12 x 2 ð35Þ where the terms M i with i ¼ 1; . . . ; 11 are explicitly given in Appendix C; they are dependent on the phase difference and the horizon. As explained in Appendix C, for long wavelength (or small wave number, j), the resulting wave dispersion relation is the same as that of the CT (Soedel, 2004). Thus, both theories give the same relationship for long wavelength. For specified values of E ¼ 200 GPa, q ¼ 7850 kg=m3 , 2

h ¼ 5  108 m, k ¼ 5=ð6  mÞ,

m ¼ 1=3, and finite horizon size,

8

d ¼ 10 m, the evaluation of the determinant without any simplifications leads to the variation of the wave frequency, x, as a function of the wave number, j. Fig. 6 shows comparisons of

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

7

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Fig. 6. Comparison of wave speed dispersions: (a) lowest flexural mode x1 , (b) thickness-shear mode x2 , and (c) thickness-twist mode x3 .

pffiffiffiffiffi 2 nondimensionalized phase speed (t=ðts k Þ) as a function of wave number j=ð2p=dÞ for the first three modes; lowest flexural mode x1 , thickness-shear mode x2 , and thickness-twist mode x3 . As observed in Fig. 6a, both the classical and PD theories estimate zero speed in the limit as wave number, j, approaches zero, whereas the classical theory phase speed nearly approaches the Rayleigh surface wave speed, t=ts ¼ 0:9274, for Poisson’s ratio of m ¼ 0:30 as wave number, j, increases (Stephen, 1997). In Fig. 6b and c, both theories estimate comparable results for the long wavelengths (or relatively small j). However, PD theory captures the feature of real materials that phase velocity decreases as the wave number increases. Also comparisons of wave frequency dispersions for increasing wave number are shown in Fig. 7. As a characteristic of real materials, dispersion curves of the peridynamic theory for all modes level off as the wave number increases and exceeds a value of 2p=d (Silling, 2000). Thus, PD theory captures the experimentally observed feature of real materials, which are always dispersive as a result of long-range forces. 6. Critical curvature and critical angle In order to include failure in the material response, the response functions in the governing equations for the plate can be modified through a history-dependent scalar value function, HðxðjÞ  xðkÞ ; tÞ as

^f ðkÞðjÞ ¼ cs HðxðjÞ  xðkÞ ; tÞðuðkÞðjÞ Þ

ð36aÞ

HðxðjÞ  xðkÞ ;tÞ   (   0  jðkÞðjÞ ðxðjÞ  xðkÞ ; t0 Þ < jc and u 1 if ðkÞðjÞ ðxðjÞ  xðkÞ ;t Þ < uc ¼ 0 otherwise ð37Þ

Critical curvature and angle values can be expressed in terms of the critical energy release rate of the material. In order to find these relationships, the total strain energy required to remove all of the interactions across a newly created crack surface, A, shown in Fig. 8, must be determined and equated to the corresponding critical energy release rate value. The total bending strain energy required to remove all of the interactions across the new crack surface A is 

W cb ¼

J Kþ X X 1 k¼1 j¼1

2

    cb ðjc Þ2 xðj Þ  xðkþ Þ  V ðkþ Þ V ðj Þ

ð36bÞ

ð38Þ

The total bending strain energy, W cb , can be equated to the mode-I critical energy release rate, GIc , in order to determine the value of the bending critical curvature as

GIc ¼

1 c ð 2 b

jc Þ2

 PK þ PJ     þ þ k¼1 j¼1 xðj Þ  xðk Þ  V ðk Þ V ðj Þ A

ð39Þ

Based on the expression derived by Silling and Askari (2005) and Madenci and Oterkus (2014) for the critical energy release rate, it is evident that

 PK þ PJ     þ þ k¼1 j¼1 xðj Þ  xðk Þ  V ðk Þ V ðj Þ

and

~f ðkÞðjÞ ¼ cb HðxðjÞ  xðkÞ ; tÞðjðkÞðjÞ Þ

It is defined as

A

¼

hd4 2

ð40Þ

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

8

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Fig. 7. Comparison of wave frequency dispersions: (a) lowest flexural mode x1 , (b) thickness-shear mode x2 , and (c) thickness-twist mode, x3 .

GIIIc ¼

1 cð 2 s

uc Þ2

 PK þ PJ     þ þ k¼1 j¼1 xðj Þ  xðk Þ  V ðk Þ V ðj Þ A

ð43Þ

By using the relationship given in Eq. (40), the critical shear angle can be obtained as

uc ¼

sffiffiffiffiffiffiffiffiffiffiffiffi 4GIIIc cs hd4

ð44Þ

7. Results

Fig. 8. Interaction between two material points whose line of action crosses the crack surface.

As part of the numerical results, simple static loading conditions are considered first to compare the PD predictions with the analytical solutions. A plate with a center crack under bending is considered next. In order to obtain the static solution, the adaptive dynamic relaxation technique given by Kilic and Madenci (2010) is used, and the horizon size is chosen as d ¼ 3:015Dx where Dx is the uniform grid spacing.

Finally, the critical curvature can be expressed as

jc ¼

sffiffiffiffiffiffiffiffiffiffiffiffi 4GIc

ð41Þ

cb hd4

Similarly, the total shear strain energy required to remove all of the interactions across the surface A is 

W cs ¼

J Kþ X X 1 k¼1 j¼1

2

    cs ðuc Þ2 xðj Þ  xðkþ Þ  V ðkþ Þ V ðj Þ

ð42Þ

This total shear strain energy, W cs , can be equated to the mode-III critical energy release rate, GIIIc , in order to determine the value of the critical shear angle as

7.1. Timoshenko beam subjected to pure bending and transverse loading The length of the beam is L ¼ 1 m, with a cross sectional area of A ¼ 0:1  0:1 m2 . Its Young’s modulus is specified as E ¼ 200 GPa. Only a single row of material (collocation) points are necessary to discretize the beam. The distance between material points is Dx ¼ 0:01 m. The left edge is constrained by introducing a fictitious region with a size of d. The beam is first subjected to bending moment and then transverse loading, as shown in Figs. 9 and 10. The loading is applied to a single material point at the right end ~ ¼ 3:33  109 N=m2 for bending and of the bar as a body load of b

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

9

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

As depicted in Fig. 11, the PD and analytical solutions are in excellent agreement.Under the transverse loading case, the analytical solutions for the transverse displacement and the rotation are given as

w¼

  Px P x3 Lx2  þ kGA 2EI 3

ð46aÞ

Pð2Lx  x2 Þ 2EI

ð46bÞ

and

/¼

As shown in Fig. 12, the PD and the analytical solutions also agree well with each other. This verifies that the PD equations of motion accurately captures the deformation behavior of a Timoshenko beam. 7.2. Mindlin plate subjected to pure bending and transverse force loading

Fig. 9. (a) Timoshenko beam subjected to pure bending and (b) its discretization.

As shown in Figs. 13 and 14, the length and width of the plate is L ¼ W ¼ 1 m with a thickness of h ¼ 0:1 m. The Young’s modulus of the plate is specified as E ¼ 200 GPa. Only a single row of material (collocation) points in the thickness direction is necessary to discretize the domain. The distance between material points is Dx ¼ 0:01 m. The left edge is constrained by introducing a fictitious region with a size of 3Dx. The loading is applied to a single row of material points at the right end of the plate as a resultant body load ~x ¼ 3:33  108 N=m for bending and b ^ ¼ 5  108 N=m2 for the of b transverse loading. The peridynamic solutions of the transverse displacement and the x-direction rotation for bending moment and transverse loading cases are compared with finite element (FE) solutions by using a shell element, which is suitable for thick shell structures, available in commercial software, ANSYS. As depicted in Figs. 15 and 16, the PD and the FE solutions agree well with each other. This verifies that the PD equation of motion given in Eqs. (30a-c) can accurately capture the deformation behavior of a Mindlin plate. 7.3. Mindlin plate with central crack Crack growth in a square plate with an initial central crack aligned with the x-axis, as shown in Fig. 17, is analyzed. The length and width of the plate are L ¼ W ¼ 1 m with a thickness of h ¼ 0:1 m. Plate thickness to crack length ratio is h=2a ¼ 0:5, which has the properties of a thick plate, where 2a is the initial crack length. The Young’s modulus of the plate is specified as E ¼ 3:227 GPa and the shear modulus is G ¼ 1:21 GPa. Only a single row of material (collocation) points in the thickness direction is necessary to discretize the domain. The distance between mate-

Fig. 10. (a) Timoshenko beam subjected to transverse loading and (b) its discretization.

^ ¼ 5  109 N=m3 for the transverse loading, corresponding to an b applied moment of M ¼ 3:33  105 Nm and a transverse load of P ¼ 5  105 N, respectively. The analytical solutions for the transverse displacement and the rotation due to pure bending are given as

Mx2 w¼ 2EI

ð45aÞ

and

/¼

Mx EI

ð45bÞ

rial points is Dx ¼ 2  103 m. The horizon size is chosen as d ¼ 3:015Dx. The material is chosen as polymethyl-methacrylate (PMMA), which shows a brittle fracture behavior. Mode-I fracture toughness pffiffiffiffiffi of this material is given as 1:33 MPa m (Ayatollahi and Aliha, pffiffiffiffiffi 2009) and Mode-III fracture toughness is given as 7:684 MPa m (Farshad and Flueler, 1998). The critical energy release rates of mode-I and mode-III can be found from

GIc ¼

K 2Ic E

and GIIIc ¼

K 2IIIc 2G

ð47Þ

In order to show simple mode-I crack growth, a bending moment loading is applied through a single row of material points at the horizontal boundary regions of the plate. Small increments ~y ¼ 250 N=m are induced in order of resultant body loading of Db

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

10

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Fig. 11. Variation of (a) rotation and (b) transverse displacement along a Timoshenko beam subjected to pure bending loading.

Fig. 12. Variation of (a) rotation and (b) transverse displacement along a Timoshenko beam subjected to transverse force loading.

to achieve stable crack growth. Under the applied uniform bending, the crack starts to grow at the end of nearly 66,000 time steps, and as expected, it propagates towards the edges of the plate, as shown in Fig. 18(a–d).

Substituting from Eqs. (A.1a,b) into (23a,b) and performing the algebraic manipulations lead to

€ ðkÞ ¼ cs qw

1   X ^ wðkÞ;xx  /ðkÞ;x nðjÞðkÞ V ðjÞ þ b ðkÞ

ðA:2aÞ

j¼1

8. Final remarks

and

This study presented the PD equations of motion for a Timoshenko beam and Mindlin plate. PD dispersion relationships were also obtained and it was observed that they are similar to the ones observed in real-materials, which cannot be predicted by using classical theory. After establishing the validity of the PD predictions, the expressions for critical curvature and shear angle values in terms of mode-I and mode-III critical energy release rates of the material were also utilized to predict crack growth in a plate under pure bending.

qI € A

/ðkÞ ¼ cb

1 X /ðkÞ;xx nðjÞðkÞ V ðjÞ j¼1

þ cs

 1  X 1 ~ wðkÞ;x  /ðkÞ  /ðkÞ;xx n2ðjÞðkÞ nðjÞðkÞ V ðjÞ þ b ðkÞ ðA:2bÞ 4 j¼1

Expressing the infinitesimal incremental volume, V ðjÞ ¼ ADnðjÞðkÞ , with DnðjÞðkÞ representing the spacing between two consecutive material points, and converting the summation to an integration as DnðjÞðkÞ approaches zero, i.e. DnðjÞðkÞ ! dn, yield

Appendix A. PD material parameters

€ ¼ cs qw

Z

d



^ w;xx  /;x nAdn þ b

ðA:3aÞ

0

A.1. Timoshenko beam and As the horizon size approaches zero, PD equations must recover their classical counterparts. Therefore, the out-of-plane deflection and transverse shear angle at material point j can be expressed, using Taylor series expansion and disregarding higher order terms, as



qI € A

/ ¼ cb

ðA:1a; bÞ

d

/;xx n Adn þ cs

Z

0

0

d

  1 ~ w;x  /  /;xx n2 nAdn þ b 4 ðA:3bÞ

Performing the integrations results in

1 wðjÞ ¼ wðkÞ þ wðkÞ;x nðjÞðkÞ sgn xðjÞ  xðkÞ þ wðkÞ;xx n2ðjÞðkÞ 2

1 /ðjÞ ¼ /ðkÞ þ /ðkÞ;x nðjÞðkÞ sgn xðjÞ  xðkÞ þ /ðkÞ;xx n2ðjÞðkÞ 2

Z

€ ¼ Acs qw

d2 ^ ðw;xx  /;x Þ þ b 2

ðA:4aÞ

and

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

11

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Fig. 13. (a) Mindlin plate subjected to pure bending loading and (b) its discretization. 2

qI € A

/¼A

d 2

2

cb  cs

!

2

d d ~ /;xx þ A cs ðw;x  /Þ þ b 8 2

ðA:4bÞ

Note that these PD equations of motion have the same form as those of the classical Timoshenko beam equations

^ € ¼ k G ðw;xx  /;x Þ þ b qw

Fig. 14. (a) Mindlin plate subjected to transverse force loading and (b) its discretization.

/aðjÞ ¼ /aðkÞ þ /aðkÞ;x nðjÞðkÞ cos h þ /aðkÞ;y nðjÞðkÞ sin h 1 þ /aðkÞ;xx n2ðjÞðkÞ cos2 h þ /aðkÞ;xy n2ðjÞðkÞ cos h sin h 2 1 2 þ /aðkÞ;yy n2ðjÞðkÞ sin h 2

ðA:5aÞ

and

qI € EI ~ /¼ /;xx þ k G ðw;x  /Þ þ b A

A

ðA:5bÞ

where the shear correction factor, k, can be taken as 5=6 for rectangular cross sections. Finally, equating the coefficients of the independent variables in the PD equations of motion to those of the classical equations yields the relationships between the PD material constants, cs and cb , and the shear and Young’s moduli, G and E, as

cs ¼

2k G Ad2

and cb ¼

2E I d2 A2

þ

1 kG 4 A

ðA:6a; bÞ

Substituting from Eqs. (A.7a,b) into (30a-c) and performing the algebraic manipulations lead to

^ p d3 € ¼ cs h qhw w;xx þ w;yy  /x;x  /y;y þ b 2

q

3

p 2

cs h

d3 ~x ðw;x þ /x Þ þ b 3

ðA:8bÞ

and

As the horizon size approaches zero, PD equations must recover their classical counterparts. Therefore, out-of-plane deflection and rotations at material point j can be expressed, using Taylor series expansion and disregarding higher order terms, as

! 3

h € p 3 2 d2 3/y;yy þ /y;xx þ 2/x;xy q /y ¼ hd cb  cs 3 12 16 10 

wðjÞ ¼ wðkÞ þ wðkÞ;x nðjÞðkÞ cos h þ wðkÞ;y nðjÞðkÞ sin h 1 þ wðkÞ;xx n2ðjÞðkÞ cos2 h þ wðkÞ;xy n2ðjÞðkÞ cos h sin h 2 1 2 þ wðkÞ;yy n2ðjÞðkÞ sin h 2

ðA:8aÞ

! 3

h € p 3 2 d2 /x ¼ hd cb  cs 3/x;xx þ /x;yy þ 2/y;xy 3 12 16 10 

A.2. Mindlin plate

ðA:7bÞ

p 2

cs h

~ d3 w;y þ /y þ b y 3

ðA:8cÞ

The classical counterpart of the PD equations of motion for a Mindlin plate can be written as

ðA:7aÞ



^ € ¼ hk2 G w;xx þ w;yy  /x;x  /y;y þ b qhw

ðA:9aÞ

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

12

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Fig. 15. Variation of (a) rotation and (b) transverse displacement along a Mindlin plate subjected to pure bending loading.

Fig. 16. Variation of (a) rotation and (b) transverse displacement along a Mindlin plate subjected to transverse force loading.

  D D /x ¼ D/x;xx þ ð1  mÞ/x;yy þ ð1 þ mÞ/y;xy 2 2 12 2 ~  k hGð/  w;x Þ þ bx

qh3 €

x

and

ðA:9bÞ

  qh € D D /y ¼ D/y;yy þ ð1  mÞ/y;xx þ ð1 þ mÞ/x;xy 2 2 12 2 ~  k hGð/  w;y Þ þ by 3

y

SED due to bending in the classical theory can be expressed as

ðA:9cÞ

2

where k is a correction coefficient that is introduced to account for the fact that the shear stresses are not constant over the thickness. The parameters D and G are the flexural rigidity and shear modulus, respectively, which are defined as 3

Eh D¼ 12ð1  m2 Þ

E and G ¼ 2ð1 þ mÞ

ðA:10a; bÞ

The correction factor can be chosen based on the frequency of the lowest thickness shear mode as 2

p2

ðA:11Þ

12

Finally, equating equations of the peridynamic and classical theories reveals the relationships between the peridynamic material constants, cs and cb , and Young’s modulus, E, as well as constraint on the value of the Poisson’s ratio as

cs ¼

Appendix B. Surface corrections B.1. Timoshenko beam

and

k ¼

m ¼ 1=3.

9E

k 4pd3

2

2

and cb ¼

E 3h 27 2 þ k pd 4d2 80

! ðA:12a; bÞ

W CM b ðxðkÞ Þ ¼

  1 EI j2 2 A

ðB:1Þ

where j ¼ /;x . Its counterpart in PD theory can be expressed in discretized form as

W PD b ðxðkÞ Þ ¼

1 1X cb ðjðkÞðjÞ Þ2 nðjÞðkÞ V ðjÞ 2 j¼1 2

ðB:2Þ

Considering pure bending loading, the PD SED becomes

W PD b ðxðkÞ Þ ¼

1 1X cb ðjÞ2 nðjÞðkÞ V ðjÞ 2 j¼1 2

ðB:3Þ

Hence, the surface correction factor for pure bending can be defined as

g b ðxðkÞ Þ ¼

W CM b ðxðkÞ Þ W PD b ðxðkÞ Þ

ðB:4Þ

On the other hand, SED due to transverse shear deformation in the classical theory can be written as

W CM s ðxðkÞ Þ ¼

1 GkðuÞ2 2

ðB:5Þ

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

13

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

B.2. Mindlin plate SED due to bending in the classical theory can be expressed as

W CM b ðxðkÞ Þ ¼

  1 1m 2 D j2x þ 2mjx jy þ j2y þ jxy 2 2

ðB:10Þ

where jx ¼ /x;x , jy ¼ /y;y , and jxy ¼ /x;y þ /y;x . Its counterpart in the PD theory can be expressed in discretized form as

W PD b ðxðkÞ Þ ¼

1 1X cb ðjðkÞðjÞ Þ2 nðjÞðkÞ V ðjÞ 2 j¼1 2

ðB:11Þ

Considering pure bending loading along the x- and y-directions, i.e.

Mx ¼ M y –0;

M xy ¼ 0

ðB:12Þ

results in curvatures

jx ¼ jy –0; jxy ¼ 0

ðB:13Þ

The classical and PD SED for this loading condition become

W CM b ðxðkÞ Þ ¼

 1  2 D jx þ 2mjx jy þ j2y 2

ðB:14aÞ

1 2 1X cb  jx cos2 ðhÞ þ jy sin2 ðhÞ nðjÞðkÞ V ðjÞ 2 j¼1 2

ðB:14bÞ

and

W PD b ðxðkÞ Þ ¼

Hence, the surface correction factor for pure bending can be defined as

g b ðxðkÞ Þ ¼ Fig. 17. (a) Mindlin plate with a central crack subjected to pure bending loading and (b) its discretization.

where u ¼ w;x  /. Its counterpart in the PD theory can be expressed in discretized form as

W PD s ðxðkÞ Þ ¼

1 1X cs

2

j¼1

2

ðuðkÞðjÞ Þ2 nðjÞðkÞ V ðjÞ

ðB:6Þ

ðB:7Þ

g s ðxðkÞ Þ ¼

W PD s ðxðkÞ Þ

g b ðxðkÞ Þ þ g b ðxðjÞ Þ 2

ðB:8Þ

ðB:9aÞ

and

g ðxðkÞ Þ þ g s ðxðjÞ Þ gs ¼ s 2

On the other hand, SED due to transverse shear deformation in the classical theory can be written as

W CM s ðxðkÞ Þ ¼

1 2 Gk hððux Þ2 þ ðuy Þ2 Þ 2

ðB:16Þ

where ux ¼ w;x  /x and uy ¼ w;y  /y . Its counterpart in the PD theory can be expressed in discretized form as 1 1X cs ðu Þ2 n V ðjÞ 2 j¼1 2 ðkÞðjÞ ðjÞðkÞ

ðB:17Þ

Applying pure shear loading along the x-direction, i.e.

Q x –0 and Q y ¼ 0

ðB:18Þ

ux –0 and uy ¼ 0

ðB:19Þ

Classical and PD SED for this loading condition become

The surface correction factors for material point j and material point k can have different values; therefore, the surface correction factors for an interaction between material point j and material point k can be taken as their average

gb ¼

ðB:15Þ

results in shear angles

Similar to pure bending, the correction factor for pure shear loading can be obtained as

W CM s ðxðkÞ Þ

W PD b ðxðkÞ Þ

W PD s ðxðkÞ Þ ¼

Considering pure shear loading along the x-direction, PD SED becomes 1 1X cs W PD ðuÞ2 nðjÞðkÞ V ðjÞ s ðxðkÞ Þ ¼ 2 j¼1 2

W CM b ðxðkÞ Þ

W CM s ðxðkÞ Þ ¼

ðB:20aÞ

1 1X cs ðu cosðhÞÞ2 nðjÞðkÞ V ðjÞ 2 j¼1 2 x

ðB:20bÞ

and

W PD s ðxðkÞ Þ ¼

Similar results can also be obtained for pure shear loading along the y-direction as

W CM s ðxðkÞ Þ ¼ ðB:9bÞ

1 2 Gk hðux Þ2 2

1 2 Gk hðuy Þ2 2

ðB:21aÞ

and

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

14

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Fig. 18. Crack evolution at (a) 66,000th, (b) 67,000th, (c) 68,000th, and (d) 69,000th time step.

W PD s ðxðkÞ Þ ¼

1 1X cs 2 ðu sinðhÞÞ nðjÞðkÞ V ðjÞ 2 j¼1 2 y

ðB:21bÞ

g s ðxðkÞ Þ ¼

ðB:22Þ

The surface correction factors for material point j and material point k can have different values; therefore, the surface correction factors for an interaction between material point j and material point k can be taken as their average

gb ¼

g b ðxðkÞ Þ þ g b ðxðjÞ Þ 2

ðB:23aÞ

and

gs ¼

g s ðxðkÞ Þ þ g s ðxðjÞ Þ 2

Z

d

ð1 þ cosðjnÞÞndn

ðC:1aÞ

ð1  cosðjnÞÞ dn n

ðC:1bÞ

i sinðjnÞdn

ðC:1cÞ

0

Similar to pure bending, the correction factor for pure shear loading can be obtained as

W CM s ðxðkÞ Þ W PD s ðxðkÞ Þ

B1 ¼

B2 ¼

Z

d 0

and

B3 ¼

Z

d 0

As suggested by Silling (2000), in the limit of long wavelength (or small j), these integrals can be analytically evaluated by considering the first three terms of the Taylor series expansion of the cosine and sine functions

cosðjnÞ ¼ 1  ðB:23bÞ

ðjnÞ2 ðjnÞ4 þ   2! 4!

ðC:2aÞ

and

ðjnÞ3 ðjnÞ5 þ   3! 5!

Appendix C. Dispersion relations

sinðjnÞ ¼ jn 

C.1. Timoshenko beam

With the evaluation of these integrals and considering the PD material parameters, the determinant from Eq. (33) can be expressed as

The terms appearing in Eq. (33) are defined as

ðC:2bÞ

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

X4pd  X2pd 2ð1 þ mÞf2 þ kf2 þ mk  2

ð1 þ mÞf4 kf4 mkf4  þ 12t 2 24t 2 576t 4

!

Xccm ¼ Xpd

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u !2 u u 2 2 u u 2 f k km f k km ¼ tf ð1 þ mÞ þ þ  t f2 ð1 þ mÞ þ þ  f4 ð2kð1 þ mÞÞ 2 2 2 2

2

kð1 þ mÞf8 mk f8 kð1 þ mÞf6 31mk f6 þ   288t4 2560t6 6t 2 2880t 4 2 10 2 4 mk f mk f þ 2kð1 þ mÞf4  þ ¼0 129; 600t8 8t 2

15

þ

ðC:7Þ

ðC:3Þ

in which m is the Poisson’s ratio. The nondimensional wave frequency and wave number, Xpd and f, respectively, are defined as

and the second mode can be expressed as

Xccm ¼ Xpd

ðC:4a; bÞ

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u !2 u u 2 2 u u f k km t 2 f k km f ð1 þ mÞ þ ¼ tf2 ð1 þ mÞ þ þ þ þ  f4 ð2kð1 þ mÞÞ 2 2 2 2

The non-dimensional geometric parameters, m and t, are defined as

ðC:8Þ

Xpd ¼

m¼

xh and f ¼ jh ts

Ah I

2

and t ¼

h d

ðC:5a; bÞ

Disregarding the higher order terms of the horizon simplifies the wave dispersion relation for long wavelength limit (or small f) as



X4pd  X2pd 2ð1 þ mÞf2 þ kf2 þ mk þ 2kð1 þ mÞf4 ¼ 0

ðC:6Þ

As expected, this equation recovers the nondimensional wave dispersion relation in the classical theory (CT) (Amirkulova, 2011) for a long wavelength limit. It leads to four different values for wave frequency, which represent two waves traveling to the right and two waves traveling to the left of the beam. Therefore, there are two distinct modes that can propagate in a Timoshenko beam (Reis, 1978). The first mode yields zero frequency (X ¼ 0) when the wave number is equal to zero (f ¼ 0)

The wave dispersion relations for long wavelength limit are shown in Figs. C.1 and C.2 for the specified values of E ¼ 200 GPa,

q ¼ 7850 kg=m3 , k ¼ 5=6, h ¼ 107 m, m ¼ 0:3, and horizon size, d ¼ 108 m. Note that the wave number of the second mode is real for X P 3:16. This indicates that the first mode is a propagating mode for any wave frequency, while the second mode only propagates for wave frequencies X P 3:16 and is exponentially attenuated for X < 3:16, as discussed by Reis (1978). Therefore, X ¼ 3:16 is the cut-off frequency for the second mode. pffiffiffi The variations of nondimensional phase speed, t ¼ X=ðf kÞ and pffiffiffiffiffiffiffiffiffi t ¼ X=ðf E=GÞ, as a function of the wave number are shown in Fig. C.3 for the first and second modes. pffiffiffi pffiffiffiffiffiffiffiffiffi The phase speed, t ¼ X=ðf kÞ or t ¼ X=ðf E=GÞ, converges to unity as the wave number, f, increases. Hence, the phase speed

Fig. C.1. Classical wave frequency dispersion for the first and second modes (combine these plots).

Fig. C.2. Classical wave speed dispersion for the first and second modes.

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

16

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

Fig. C.3. Classical phase speed dispersions and wave frequency dispersions.

pffiffiffi for the first mode, t ¼ X=ðf kÞ, converges to the shear wave speed, ts , of the CT, since

t¼

x X t pffiffiffi ¼ pffiffiffi and j f k ts k

Moreover, the pffiffiffiffiffiffiffiffiffi

ðC:9Þ

phase

speed

for

the

second

mode,

t ¼ X=ðf E=GÞ, converges to the compressional wave speed, tc , of the CT, since

t

X

qffiffiffi ¼ f GE tc

ðC:10Þ

M11 ¼

Z 2p Z 0

0

d

eijn cosðhÞ þ 1 2 nsin ðhÞdV 2

in which dV denotes the infinitesimal volume of a material point, which can be written in cylindrical coordinates as dV ¼ hndndh. Evaluation of these integrals yields Bessel functions of the first kind, J0(jd) and J1(jd), and Struve functions, H0(jd) and H1(jd). As indicated by Silling (2000), in the limit of a long wavelength (C) or for a very small wave number (j ! 0), the integrals in Eqs. (C.11a-k) can be simplified by using the first three terms of the Taylor series expansion for cosine and sine functions

cosðjn cosðhÞÞ ¼ 1  C.2. Mindlin plate

ðjn cosðhÞÞ2 ðjn cosðhÞÞ4 þ   2! 4!

M1 ¼

0

M2 ¼

0

Z 2p Z

Z 2p Z

Z 2p Z

Z 2p Z 0

M6 ¼

Z 2p Z

eijn cosðhÞ þ 1 sinðhÞdV 2

ðC:11cÞ

eijn cosðhÞ  1 cos2 ðhÞdV n

ðC:11dÞ

e



d

0

Z 2p Z

d

0

Z 2p Z 0

M 10 ¼

ðC:11bÞ

n d

Z 2p Z

0

M9 ¼

eijn cosðhÞ þ 1 cosðhÞdV 2

1

0

Z 2p Z 0

0

d

sinðhÞ cosðhÞdV

ðC:11eÞ

eijn cosðhÞ  1 cosðhÞdV

ðC:11fÞ

eijn cosðhÞ þ 1 n cos2 ðhÞdV 2

ðC:11gÞ

eijn cosðhÞ þ 1 n sinðhÞ cosðhÞdV 2

ðC:11hÞ



Substituting the relationships given in Eqs. (C.12a,b) into integrations given in Eqs. (C.11a-k) and solving the determinant equation given in Eq. (35), while ignoring higher order terms of horizon size, yield the dispersion relationship for a long wavelength limit in the PD theory   2 2   k Ghj2 þ qhx2 k Ghij 0     3 2 2 qh 2 2 ¼0  k Ghij Dj  k Gh þ 12 x 0     3 2  0 0  D ð1  mÞj2  k Gh þ qh x2  12

ðC:13Þ

eijn cosðhÞ  1 2 sin ðhÞdV n

d

ðjn cosðhÞÞ3 ðjn cosðhÞÞ5 þ   3! 5! ðC:12bÞ

2

0

0

M8 ¼

ðC:11aÞ

ijn cosðhÞ

d

sinðjn cosðhÞÞ ¼ jn cosðhÞ 

eijn cosðhÞ  1 dV n

0

0

M7 ¼

d

0

0

M5 ¼

d

0

0

M4 ¼

d

0

0

M3 ¼

d

ðC:12aÞ

and

The terms appearing in Eq. (35) are defined as

Z 2p Z

ðC:11kÞ

eijn cosðhÞ  1 sinðhÞdV

with the constraint on Poisson’s ratio (m ¼ 1=3). This relationship is equivalent to the dispersion relationship obtained from the classical theory (CT). Moreover, roots of Eq. (C.13) correspond to three different natural frequencies. Soedel (2004) explained that the lowest of these frequencies is the one that the transverse deflection mode dominates and other two are considered as shear modes. Shown in Fig. C.3 are the nondimensionalized phase speed (t=ts ) dispersion relationships with the change of wave number (h=C) for three different wave modes for the long wavelength limit 2

ðC:11iÞ

ðC:11jÞ

while considering the following properties of a plate k ¼ 5=ð6  mÞ and m ¼ 1=3. Variations similar to those in Fig. C.3 and comparisons with Rayleigh–Lamb waves, which have the property of waves in a plate with infinite extent, can be seen in Stephen (1997). Also, in Stephen (1997), x1 is named as the lowest flexural mode, x2 as the thickness-shear mode, and x3 as the thickness-twist mode.

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040

C. Diyaroglu et al. / International Journal of Solids and Structures xxx (2015) xxx–xxx

In Fig. C.3, when the wave number (h=C) reaches the value of two, which means the wave length is now comparable with the thickness of a plate, all the phase speeds become flat. Fig. C.3 shows the nondimensionalized wave frequency (ðt=ts Þðh=CÞ) dispersions with the change of wave number (h=C) for three wave modes of the CT or for the PD theory in the long wavelength limit. References Agwai, A., Guven, I., Madenci, E., 2011. Predicting crack propagation with peridynamics: a comparative study. Int. J. Fract. 171, 65–78. Amirkulova, F.A., 2011. Dispersion Relations for Elastic Waves in Plates and Rods (M.Sc Thesis), Rutgers, The State University of New Jersey. Ayatollahi, M.R., Aliha, M.R.M., 2009. Analysis of a new specimen for mixed mode fracture tests on brittle materials. Eng. Fract. Mech. 76, 1563–1573. Eringen, A.C., 1972. Linear theory of nonlocal elasticity and dispersion of plane waves. Int. J. Eng. Sci. 10 (5), 425–435. Farshad, M., Flueler, P., 1998. Investigation of mode III fracture toughness using an anti-clastic plate bending method. Eng. Fract. Mech. 60, 597–603.

17

Kilic, B., Madenci, E., 2010. An adaptive dynamic relaxation method for quasi-static simulations by using peridynamic theory. Theor. Appl. Fract. Mech. 53, 194– 204. Madenci, E., Oterkus, E., 2014. Peridynamic Theory and Its Applications. Springer, New York. O’Grady, J., Foster, J., 2014a. Peridynamic beams: a non-ordinary, state-based model. Int. J. Solids Struct. 51, 3177–3183. O’Grady, J., Foster, J., 2014b. Peridynamic plates and flat shells: a non-ordinary, state-based model. Int. J. Solids Struct. 51, 4572–4579. Reis, M. de A. e S., 1978. Wave Propagation in Elastic Beams and Rods (Ph.D. Dissertation), Massachusetts Institute of Technology. Silling, S.A., 2000. Reformulation of elasticity theory for discontinuities and longrange forces. J. Mech. Phys. Solids 48, 175–209. Silling, S.A., Askari, E., 2005. A meshfree method based on the peridynamic model of solid mechanics. Comput. Struct. 83, 1526–1535. Soedel, W., 2004. Vibrations of Shells and Plates, third ed. Marcel Dekker, New York. Stephen, N.G., 1997. Mindlin plate theory: best shear coefficient and higher spectra validity. J. Sound Vib. 202, 539–553. Taylor, M., Steigmann, D.J., 2013. A two-dimensional peridynamic model for thin plates. Math. Mech. Solids, 1–13. Weckner, O., Silling, S.A., 2011. Determination of nonlocal constitutive equations from phonon dispersion relations. Int. J. Multiscale Comput. Eng. 9 (6), 623–634.

Please cite this article in press as: Diyaroglu, C., et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int. J. Solids Struct. (2015), http://dx.doi.org/10.1016/j.ijsolstr.2015.04.040