Peridynamic simulation to impacting damage in composite laminate

Composite Structures 138 (2016) 335–341

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Peridynamic simulation to impacting damage in composite laminate Chaoyang Sun, Zaixing Huang ⇑ State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, PR China

a r t i c l e

i n f o

Article history: Available online 12 December 2015 Keywords: Peridynamics PD rate-dependent constitutive model Composite laminates Impact velocity Damage

a b s t r a c t Based on combination of the prototype microelastic brittle (PMB) model and the viscoelastic Kelvin–Voigt model, we propose a peridynamic (PD) rate-dependent constitutive equation and a new interlayer bond describing interlayer interactions of fiber reinforced composite laminate. They are used to simulate total process from low-velocity collision to high-velocity penetration appearing in composite laminate. The influences of impacting velocity on the damage pattern and damaged area of composite laminate are analyzed and discussed. The results reveal that impacting damage in composite laminate is caused by the stretching strain. With the impacting velocity increasing gradually, the damage firstly occurs at the fixed edges, subsequently on the bottom and finally on the top surface of composite laminate. Meanwhile, owing to competition between damage and damping, the damaged area of composite laminate firstly rises with an increase in impacting velocity and reaches a peak, and then decreases rapidly. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Fiber reinforced composite laminate is widely used in the aircraft industry due to its high strength, modulus and fatigue resistance. However, severe damage may occur in composite laminates due to impact loadings during their manufacture and utilization. There are various impacting forms, such as lowvelocity free falling and high-speed striking, which can result in entirely different consequences. Cantwell and Morton conducted a series of low and high velocity impact tests on carbon fiber reinforced laminates [1]. The experimental observations revealed that under low-velocity collision, the response of target is determined by its structural geometry; and under high-velocity impact, a localized form of failure is generated in the target. The experiment performed by Choi investigated the low-velocity impact damage in T300/976 graphite/epoxy laminated composites [2]. Morye developed a simple energy absorption model for the ballistic damage of polymer composites and found a good correlation between the damage radius measured by high speed photography and theoretical prediction based on the transverse wave speed [3]. The finite element method (FEM) has been used to simulate failure in composite laminates subjected to impact [4–9]. However, FEM requires that the displacement field of body should be continuously differential, but this requirement is in contradiction with inherent discontinuity existing in fracture and damage. In order ⇑ Corresponding author. E-mail address: [email protected] (Z. Huang). http://dx.doi.org/10.1016/j.compstruct.2015.12.001 0263-8223/Ó 2015 Elsevier Ltd. All rights reserved.

to solve some problems with discontinuity, Silling and his colleagues proposed and developed the so-called peridynamic theory [10]. Unlike FEM, the PD theory reformulates the basic equations of continuum in the form of integration and allows the displacement field to be discontinuous. Thus the PD theory has an advantage in dealing with the discontinuous problems, such as the crack initiation and propagation and the damage issues. Recently, the PD theory has been widely applied to analyze damage in composite laminates. Bobaru put forth a peridynamic 3D model and used it to study the complex fracture problem of carbon nanotube-reinforced polymer composites [11]. A prediction for the delamination and matrix damage of the graphite-epoxy laminates was made in detail based on the PD theory by Xu et al. [12]. Kilic employed the PD theory to predict the damage of the center-cracked composite laminates with different fiber orientations, and the results agreed with the experimental observations in the literature [13]. Oterkus combined the PD theory and FEM to calculate failure in a stiffened composite curved panel with a central slot. The results demonstrate that the PD approach may be used to assess the durability of complex composite structures [14]. Wang and Huang advanced tensor-like pairwise force function in PD theory and performed simulations to the isotropic circular plate and the anistropic lamina under the impact loading condition [15]. The studies above partly revealed the mechanism of the impact damage of composites under low or high velocity impact from different angle of view. However, a simulation to the impacting damage change from low-velocity to high-velocity

336

C. Sun, Z. Huang / Composite Structures 138 (2016) 335–341

in composite laminate still awaits investigation. This simulation relies on an appropriate constitutive model. In this paper, a PD rate-dependent constitutive equation characterizing the pairwise force is proposed based on the PMB model [16] and the viscoelastic Kelvin–Voigt model. The interlayer bond is introduced to establish the 3D models for the fiber reinforced composite laminates. The force constant of interlayer bond is derived based on the PD theory and mechanics of composite material. Using the improved LAMMPS molecular dynamic code [17], we simulate the damage of composite laminates subjected to impact with various velocities. Some new results are given. 2. Basic formulations 2.1. Peridynamic theory The peridynamic theory reformulates the basic equations of motion so that they can be applied to the problems with inherent discontinuity. In peridynamic model of elastic theory, internal forces between material points are represented as the form of integral, rather than differential formulation. As shown in Fig. 1, the equation of motion of a material point at x in the reference configuration R at time t is written as [10]

qðxÞu€ ðx; tÞ ¼

Z

fðg; nÞdV x0 þ bðx; tÞ 8x 2 R;

Fig. 2. Initial and deformed configuration in peridynamics.

ð1Þ

H

where q is the material density of the material point at x; u is the displacement vector field of the material point at x at time t; and b the prescribed body force. H is the neighborhood of the material point at x, and it is determined by a scalar d called the horizon. Within this neighborhood, the action exerted on the material point at x by another material point at x0 is characterized by the pairwise force vector f, which is a function of the relative position vector n in the reference configuration and the relative displacement vector g at time t. As shown in Fig. 2, n and g are expressed as

n ¼ x0  x;

g ¼ uðx0 ; tÞ  uðx; tÞ;

ð2Þ

According to the bond-based peridynamic theory, the pairwise force vector f (seeing Fig. 3) can be expressed as [10]

fðg; nÞ ¼ f ðg; nÞ

gþn ; jjg þ njj

ð3Þ

where f is a scalar-valued pairwise force, and |||| is the Euclidean norm. Apparently, g + n is the current relative position vector of the material points at x and x0 .

Fig. 3. The pairwise force vector in the bond-based peridynamics theory.

According to the PMB model [18], f is written as f = cs, where c = 18K/pd4, called the micro-modulus. K is the bulk modulus and s the bond stretch, which is defined as

sðg; nÞ ¼

jjg þ njj  jjnjj ; jjnjj

ð4Þ

where ||n|| is the original bond length in the reference configuration, ||g + n|| is the deformed bond length. Therefore, the bond stretch is dimensionless. Bonds break when they are stretched beyond a limitation. The critical stretch s0 of the PMB materials is given by [16]

rffiffiffiffiffiffiffiffiffi 5G0 s0 ¼ 9Kd

ð5Þ

where G0 is the critical energy release rate of material. However, an assumption implied in Eq. (5) is that the critical stretch of a bond is completely independent of other bonds. This is an oversimplification for real materials. Thus, a new timedependent critical stretch is introduced, which read [16]

s0 ðtÞ ¼ s00  asmin ðtÞ; Fig. 1. Schematic of peridynamic interactions.

ð6Þ

where s00 and a are material-dependent constants and smin(t) is the minimum stretch in all bonds connected to the material point x.

337

C. Sun, Z. Huang / Composite Structures 138 (2016) 335–341

In order to give the critical stretch criterion [18], a scalar history-dependent function l is defined as follows

lðt; g; nÞ ¼



1 if sðt 0 ; g; nÞ < s0 0

for all 0 6 t 0 6 t;

otherwise:

uðx; tÞ ¼ 1 

H

lðt; g; nÞdV x0 R H

dV x0

:

ð8Þ

From Eq. (3), it is easy to see that the pairwise force f depends only on the stretch s. Therefore, it is inadequate using the PMB model to describe impacting responses with rate-dependent effects. In order to characterize the rate effect of impact damage, we use the Kelvin–Voigt model to give an extended PMB equation, which reads

(

gþn cðs jjgþnjj

þ v s_ Þ s < s0 ;

0

s P s0 :

ð9Þ

where v is the viscous damping coefficient, and s_ is the stretch rate. Using Eq. (4), we have

s_ ¼

g_  ðg þ nÞ : jjnjj  jjg þ njj

ð10Þ

Eq. (9)is called the peridynamic rate-dependent constitutive model, which can be regarded as a simplification of the Eq. (2.27) in [19]. 3. Numerical solution In a reference configuration, a body is uniformly discretized into nodes with a certain volume [16]. Different from the finite element, all nodes form a grid allowing deformation without restriction. After discretization, the integral in Eq. (1) can be replaced by a finite sum, that is

qu€ ni ¼

ð15Þ

4. Laminates modeling

2.2. Peridynamic rate-dependent constitutive model

fðg; nÞ ¼

lðt; g; nÞV p : PN p¼1 V p

p¼1

ð7Þ

Based on the Eq. (7), the damage at a point x can be define as [16]

R

PN

uðx; tÞ ¼ 1 

X n fðunp  uni ; xp  xi ÞV p þ bi ;

4.1. PD modeling of a lamina In an isotropic PMB model, the interactions within material are direction-independent. However, if a fiber-reinforced composite lamina is concerned, the directional dependence has to be considered. Therefore, in order to model the fiber reinforced composite lamina within which the orientation of fiber is h, two kinds of PD bonds (fiber bond and matrix bond) are necessary to be introduced. As shown in Fig. 4, the fiber bond characterizes interaction between two material points along the fiber direction, whereas the matrix bond represents interaction between two material points orthogonal to the fiber direction. The material parameter introduced by the fiber bonds is cf, and the material parameter associated with the matrix bonds is cm. On the basis of Oterkus’s works [14], we take the volume fraction of fiber into account. Thus, the bond parameters cf and cm take the form below

cf ¼ 

2jE1 ðE1  E2 Þ ; PQ E1  19 E2 q¼1 nqi V q

8ð1  jÞE1 E2 ; cm ¼   E1  19 E2 pld3

where E1, E2, G12 and t12 are the engineering elastic constants of lamina, jthe volume fraction of fiber, and l the thickness of lamina;. the symbol nqi denotes the initial distance between the point i and the point q, Vq is the volume of the material point q, and Q the total number of the material points along the fiber direction. By cf and cm, the micro-modulus of any orientation is represented as c = cf cos2 / + cm sin2 /, where / is the angle between fiber and bond. 4.2. Critical stretch of bond The lamina is strong anisotropic. So its strength is different at different direction. The strength at different direction can be determined in terms of the fiber bond and matrix bond. Fig. 5 illustrates

ð11Þ

p

where f is the pairwise force, n is the time step number, the subscripts denote the node numbers, Vp is the volume of the node p, and

uni ¼ uðxi ; tn Þ:

ð12Þ

In Eq. (11), the summation runs over all points within the neighborhood {||xp  xi|| 6 d}, and the acceleration at the left side can be expressed as an explicit central difference form:

€ ni ¼ u

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

ð13Þ

where Dt is a constant time step that satisfies the stability condition given by Silling and Askari [16]:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u 2q ; Dt < tPN c p¼1 jjxp xi jj V p

ð14Þ

where N is the total number of nodes inside the horizon of the node xi. Similarly, Eq. (8)should be written as a finite sum form:

ð16Þ

Fig. 4. Two kinds of PD bonds in a lamina.

338

C. Sun, Z. Huang / Composite Structures 138 (2016) 335–341

Fig. 7. Schematic of the interlayer response of composite laminate.

Fig. 5. Force-stretch relations of the fiber bond and the matrix bond.

the force-stretch relations of the fiber bond and the matrix bond. In Fig. 5, sf and sm are the critical stretches characterizing failure of the fiber bond and matrix bond. The critical stretches of the fiber bond and the matrix bond can be expressed as [16]

sffiffiffiffiffiffiffiffiffiffiffiffi 10Gf ; sf ¼ pc f d 5

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 10Gm ; sm ¼ pcm d5

ð17Þ

where Gf and Gm are the critical energy release rates of fiber and matrix, respectively. 4.3. PD modeling of laminates Fig. 8. Geometric model of the plate subjected to impact.

Fiber reinforced composite laminates consist of multiple plies of laminas with different fiber orientations. In order to model the interlayer interaction, we introduce a new bond and call it the interlayer bond. As shown in Fig. 6, the interlayer bond is assumed to only occur between two adjacent layers. Hence the material point of Ply(k) can’t interact directly with the points of Ply(k2) and Ply(k+2). The failure of the interlayer bond means the delamination of a composite laminate. The interlayer bond constant cin can be determined based on the normal deformation between adjacent layers. As shown in Fig. 7, a homogeneous strain f is applied on the composite laminate along the thickness direction. The material points i and j belong to the adjacent layers Ply(k) and Ply(k1), respectively. The point j is in the neighborhood of the point i. So the stretch of interlayer bond between i and j can be expressed as

ðl þ lfÞ= cosðaij  bÞ  l= cos aij ð1 þ fÞ cos aij ¼  1; l= cos aij cosðaij  bÞ

kij ¼

ð18Þ

where l is the thickness of lamina, aij is the orientation angle of interlayer bond before deformation, and b is the incremental orientation angle of interlayer bond after deformation. According to the sine theorem, it is easy to give



b ¼ arctan

f sin aij cos aij ; 1 þ f cos2 aij

ð19Þ

Substituting Eqs. (19) into (18), and then expending (18) to the first-order term of f, we have

kij ¼ f cos2 aij ;

ð20Þ

Therefore, the internal force along the thickness direction of laminate can be expressed as the sum of pairwise forces given by all interlayer bonds between Ply(k1) and Ply(k). Meanwhile, this internal force is also equal to the transverse elastic force in composite laminate. Thus,

XX cin kij cos aij V j V i ¼ Em fA; i

ð21Þ

j

where Em is the elastic modulus of the matrix material, A is the projected area of the laminate in the thickness direction. Hence,

cin ¼ P P i

j

Em A ; cos3 aij V j V i

ð22Þ

The critical stretch of the interlayer bond is denoted by sin. Similar to Eq. (17), it reads

sin ¼ Fig. 6. Peridynamic model of composite laminate.

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 10Gm

pcin d5

:

ð23Þ

339

C. Sun, Z. Huang / Composite Structures 138 (2016) 335–341 Table 1 Impact damage of laminates subjected to a 60 m/s impact. Surfaces

Isotropic plate

(0°)6 laminate

(±45°/0°)s laminate

(0°)6 laminate

(±45°/0°)s laminate

Top

Bottom

Table 2 Impact damage of laminates subjected to a 100 m/s impact. Surfaces

Isotropic plate

Top

Bottom

5. Numerical simulations Based on the peridynamic rate-dependent constitutive model, we investigate the damage of composite laminate under different impact loadings. The numerical simulations are implemented via the LAMMPS molecular dynamic code [17]. As shown in Fig. 8, an isotropic circular plate, with radius of 20 mm and thickness of 3 mm, is clamped at the edge. The plate is assumed to be made of epoxy material, whose density, the elastic modulus E, bulk modulus K and Poisson’s ratio t are 1300 kg/m3, 3.792 GPa, 2.528 GPa and 0.25, respectively. The viscous damping coefficient v takes 1  107 s. The critical energy release rate GIc = 2.37  103MPa m. In terms of Parks [18], we have c = 2.8611  1021 N/m6, a = 0.25 and s00 = 0.0186. According to Eq. (14), we take the time increment Dt = 2  107 s. The composite laminates are the same as the isotropic plate in geometry and boundary conditions. They are specified to be made of carbon-epoxy material. Each of composite laminate consists of six plies of laminas with thickness of 0.5 mm. Their material density is 1580 kg/m3, and the engineering constants E1 = 142.1 GPa, E2 = 8.73 GPa, t12 = 0.33. The fiber volume percentage j is 10%. By Eqs. (16) and (20), we have cf = 3.5810  1022 N/m6, cm = 1.194  1022 N/m6 and cin = 1.433  1022 N/m6. The critical stretches of the fiber bond, the matrix bond and the interlayer bond are sf = 0.01824, sm = 0.00912 and sin = 0.00833, respectively. The

viscous damping coefficient v takes 1  107 s and the time step 5  108 s. A rigid spherical projectile with density of 7900 kg/m3 and radius of 5 mm collides at the center of the circular plate along the thickness direction, and produces impacting effects. As a target model, the circular plate is discretized into a cubic lattice structure containing 30,150 particles. The volume of each particle V = 1.25  1010 m3. The lattice constant takes 0.5 mm and horizon 1.5 mm. By means of the LAMMPS molecule dynamic code, we simulate the damage of the isotropic circular plates, the (0°)6 laminates, and the (±45°/0°)s laminates under various impacting velocities, as shown in Table 1–3. Table 1 shows the damage in three plates under the impacting velocity of 60 m/s. In the isotropic plate, it can be observed that damage occurs at the edge and on the top and bottom surface simultaneously. The damage on the bottom surface displays a cross-star pattern. In theory, geometrical shape of damage in a circular isotropic plate should be axisymmetric due to the fact that this plate is subjected to axisymmetric load and constraint. However, since we adopt discrete cubic lattice structures in the numerical simulation, the pattern of damage exhibits a cubic symmetry. From Table 1 it is easy to see that, in the isotropic plate, there are three damaged areas: edge, top and bottom surface; in the (0°)6 laminate, there are two damaged areas: edge and bottom

340

C. Sun, Z. Huang / Composite Structures 138 (2016) 335–341

Table 3 Impact damage of laminates subjected to a 200 m/s impact. Surfaces

Isotropic plate

(0°)6 laminate

(±45°/0°)s laminate

Top

Bottom

Fig. 9. The variation of damage area with impacting velocity in (±45°/0°)s laminates.

surface; and in the (±45°/0°)s laminate, there is only a damaged area: edge. These differences show on one hand that the strength of the (±45°/0°)s laminate is greater than that of the (0°)6 laminate and the strength of the (0°)6 laminate is greater than that of the isotropic plate, and on the other hand that the damage firstly

occurs at the edge, subsequently on the bottom surface and finally on the top surface. When the impacting velocity arrives at 100 m/s, the damage in three plates is shown in Table 2. Compared Table 2 with Table 1, it is clear that the damaged area increases with a rise in the impacting velocity. Table 2 exhibits that the damaged area at center on the bottom is always greater than that on the top. More interestingly, the damage on the top distributes in a ring-shaped zone. It is just a stretched zone around the colliding contact region between projectile and plate. As thus, we see that all damage occur in the stretched section. Therefore, the impacting damage in plate is caused by the stretching strain. Figures in Table 3 illustrate the damaged state of three plates subjected to high-velocity impact of 200 m/s. Under this case, the three plates have been penetrated by the projectile. However, it is easy to find that, in this mode of penetrating damage, the damage at edge decreases compared with the damage in Table 2. This is due to fact that more energy has been dissipated by penetration and damping before reaching the edge. In order to acquire information of damage changing with impacting velocity, we perform a series of simulations for the (±45°/0°)s laminates subjected to impact with various velocities ranging from 30 m/s to 200 m/s. The results are shown in Fig. 9. It is easy to see that, in the bottom surface, the damaged areas firstly increase with the impacting velocity increasing and reach a peak near to 170 m/s, and then decrease rapidly. Similar trend also occurs in the top surface. This phenomenon is due to the rate

Table 4 Comparison of the impact damage of (±45°/0°)s laminates subjected to various velocity impact by means of the PMB model and the rate-dependent model. Models Rate-dependent

PMB

60 m/s

100 m/s

140 m/s

C. Sun, Z. Huang / Composite Structures 138 (2016) 335–341

effect of material. In the low velocity impact stage, the energy dissipated by damage is dominant. With the impacting velocity rising, the damping increases. The energy dissipated by damping exceeds the energy dissipated by damage. So the change of damaged area with impacting velocity exhibits the tendency as shown in Fig. 9. From Fig. 9, it is easy to see that when the velocity ranges from 30 m/s and 160 m/s, the total damaged area of the top surface goes beyond that of the bottom surface on the whole. However, when the impacting velocity reaches 170 m/s, the total damage on the top surface becomes smaller than that on the bottom surface. The reason still consists in the rate-dependent damping effect. When the impacting velocity exceeds 170 m/s, the impacting energy has been dissipated by damping before arriving at the edge. This causes the damage decreasing on the top surface of edge. The results given by Fig. 9 are qualitatively consistent with the experimental result of Cantwell and Morton [1]. As shown in Table 4, we compare the peridynamic ratedependent model with PMB model by means of simulation to the damage responses of the (±45°/0°)s laminates subjected to impact with different velocity. It is very clear that, under the same impacting velocity, the damage predicted by the PMB model is more severe than that predicted by the rate-dependent model. The reason to produce this difference consists in that the rate-dependent model has a better energy-absorbing ability near the impact point due to the appearance of damping. This makes that the impacting energy is partly dissipated by the damping, rather than wholly used to produce the damage. 7. Conclusions In this paper, a PD rate-dependent constitutive equation is proposed based on combination of the PMB model and the viscoelastic Kelvin–Voigt model. Meanwhile, a new interlayer bond is introduced to describe interlayer interactions in the fiber reinforced composite laminates. By means of the LAMMPS molecular dynamic code, a series of simulations are performed. Some conclusions are given as follows. (1) Peridynamics can be conveniently used to simulate total process from low-velocity collision to high-velocity penetration of composite laminate. The results show that the composite laminate has the advantage of the impact resistance compared with the isotropic plate. (2) Impacting damage in composite laminate is caused by the stretching strain. With the impacting velocity increasing gradually, the damage firstly occurs at the edge, subsequently on the bottom and finally on the top surface. (3) Owing to competition between damage and damping, the damaged area of composite laminate firstly rises with an increase in impacting velocity and reaches a peak, and then decreases rapidly. (4) Under the same impacting velocity, the damage predicted by the PMB model is more severe than that predicted by the rate-dependent model.

341

Acknowledgements The authors are grateful for the support of this work by National Natural Science Foundation of China (No. 11172130) and the Research Fund of State Key Laboratory of Mechanics and Control of Mechanical Structures (Nanjing University of Aeronautics and Astronautics) (Grant No. 0215G01). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.compstruct.2015. 12.001. References [1] Cantwell WJ, Morton J. Comparison of the low and high velocity impact response of CFRP. Composites 1989;20(6):545–51. [2] Choi HY, Downs RJ, Chang FK. A new approach toward understanding damage mechanisms and mechanics of laminated composites due to low-velocity impact: part I—experiments. J Compos Mater 1991;25(8):992–1011. [3] Morye SS, Hine PJ, Duckett RA, Carr DJ, Ward IM. Modelling of the energy absorption by polymer composites upon ballistic impact. Compos Sci Technol 2000;60(14):2631–42. [4] Hou JP, Petrinic N, Ruiz C, Hallett SR. Prediction of impact damage in composite plates. Compos Sci Technol 2000;60(2):273–81. [5] Olsson R. Analytical prediction of large mass impact damage in composite laminates. Compos Part A-Appl S 2001;32(9):1207–15. [6] Turon A, Camanho PP, Costa J, Dávila CG. A damage model for the simulation of delamination in advanced composites under variable-mode loading. Mech Mater 2006;38(11):1072–89. [7] Lim CT, Shim VPW, Ng YH. Finite-element modeling of the ballistic impact of fabric armor. Int J Impact Eng 2003;28(1):13–31. [8] Aminjikarai SB, Tabiei A. A strain-rate dependent 3-D micromechanical model for finite element simulations of plain weave composite structures. Compos Struct 2007;81(3):407–18. [9] Gu X, Xu X. Numerical simulation of damage in fiber reinforced composite laminates under high velocity impact. Acta Mater Compos Sin 2012;29 (1):150–61. [10] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 2000;48(1):175–209. [11] Bobaru F, Silling SA. Peridynamic 3D models of nanofiber networks and carbon nanotube-reinforced composites. MATERIALS PROCESSING AND DESIGN: modeling, simulation and applications-NUMIFORM 2004-proceedings of the 8th international conference on numerical methods in industrial forming processes, 712(1). AIP Publishing; 2004. p. 1565–70. [12] Xu J, Askari A, Weckner O, Silling SA. Peridynamic analysis of impact damage in composite laminates. J Aerospace Eng 2008;21(3):187–94. [13] Kilic B, Agwai A, Madenci E. Peridynamic theory for progressive damage prediction in center-cracked composite laminates. Compos Struct 2009;90 (2):141–51. [14] 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 (3):839–50. [15] Wang F, Huang Z. Peridynamic model and analysis for impact damage of composite laminates. Chin J Comput Mech 2014;31(6):709–13. [16] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Comput Struct 2005;83(17):1526–35. [17] Parks ML, Lehoucq RB, Plimpton SJ, Silling SA. Implementing peridynamics within a molecular dynamics code. Comput Phys Commun 2008;179 (11):777–83. [18] Parks ML, Plimpton SJ, Lehoucq RB, Silling SA. Peridynamics with LAMMPS: A user guide. Technical Report SAND 2008–1035, Sandia National Laboratory; January 2008. [19] Kilic B. Peridynamic theory for progressive failure prediction in homogeneous and heterogeneous materials [Ph.D. Thesis], University of Arizona; 2008.