Discrete element modeling of the microbond test of fiber reinforced composite

Computational Materials Science 49 (2010) 253–259

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Discrete element modeling of the microbond test of fiber reinforced composite Dongmin Yang a,b, Yong Sheng a,*, Jianqiao Ye a, Yuanqiang Tan b a b

School of Civil Engineering, University of Leeds, Leeds LS2 9JT, UK School of Mechanical Engineering, Xiangtan University, Hunan 411105, China

a r t i c l e

i n f o

Article history: Received 3 February 2010 Received in revised form 18 March 2010 Accepted 4 May 2010

Keywords: DEM Microbond test Fiber reinforced composite Interface debonding Matrix cracking

a b s t r a c t Discrete element method (DEM) was used to simulate the dynamic process of microbond test of fiber reinforced composites, in which the fiber and matrix were modeled as elastic and elastic–plastic materials, respectively. The interface between fiber and matrix was represented by a bilinear contact softening model. Plastic deformation and progressive cracking of the matrix were observed in the simulation with comparable similarity to the existing experimental results. The initiation and propagation of interfacial debonding were also captured by the DEM simulations, whereas it is very difficult to achieve this by other numerical methods. Vertical and inclined vises with two different vise angles were tested in the simulations. It was found that both vise geometry and vise angle had significant effect on the damages of the material, and the inclined vise was sensitive to the vise angle in terms of large variation of cutting force. The developed DEM model can also be applied to predict material damages in other more complicated fiber reinforce composite system. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Fiber reinforced composites have a number of attractive material properties characterized by their high strength and stiffness to mass ratios, damage tolerance and corrosion resistance, which make them suitable for use in many structural applications such as wind turbine, armor, naval and aerospace structures [1]. A major concern in the use of composite materials is their susceptibility to damage resulting from complicated external loading. The inherent multiphase microstructure of composites, i.e., fiber volume fraction, fiber direction and interfacial strength, results in complex failure models dominated by fractures of constituent elements [2,3]. The damage evolution and failure mechanism of composites are obviously far more complex than those of monolithic materials. They are complex processes developing simultaneously on all internal length scales, i.e., from micro to macro scale [4,5]. Single fiber composite tests, e.g. pull-out test, push-out test, and microbond test, are the widely used methods to investigate the damage of composite constituents and to study the material interactions [6–8]. Numerical modeling of microbond test, especially finite element method (FEM), has been carried out to study the interfacial strength [9–11]. The results from FEM simulation were then utilized to assess the interfacial properties, especially interface shear strength. However, matrix cracking, and its effect on interfacial debonding are very difficult to account for in FEM simulations.

* Corresponding author. Tel.: +44 (0)113 3433200; fax: +44 (0)113 3432265. E-mail address: [email protected] (Y. Sheng). 0927-0256/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2010.05.003

Normally, the propagation path of matrix cracking has to be predefined according to experimental investigations. In this paper, a discrete element method (DEM) with a commercial code of PFC2D [12] is adopted to simulate the dynamic process of microbond test of a glass fiber reinforce epoxy matrix composite system, due to its intrinsic characteristics and advantages in modeling micro scale fracturing and failure events. The objective of this study is not only to reproduce and analyze the damage process of the microbond test, but also more importantly to validate the application of DEM in fiber reinforce composite materials for further development of large scale models. In Section 2 of this paper, the basic principle of DEM will be explained, and the construction of the DEM model for composite material will be proposed in Section 3. The simulations of microbond tests are described in Section 4 and the results are analyzed in Section 5. Conclusions are drawn in Section 6.

2. Discrete element method Discrete (or distinct) element method (DEM) was firstly proposed by Cundall in 1979 [13] to study the discontinuous mechanical behavior of rock. DEM represents the material through discrete rigid elements that can interact with their neighboring ones. The rigid elements are usually small particles of disc in 2D or spherical in 3D, due to the simplicity of contact algorithm and computational advantages. Unlike the traditional finite difference method (FDM) and finite element method (FEM), which are all based on continuum mechanics, DEM can simulate the granular material with

254

D. Yang et al. / Computational Materials Science 49 (2010) 253–259

large deformation as well as the fracture behavior of solid material from continuum to discontinuum, especially brittle materials like concrete [14,15], rock [16–18] and ceramic [19–21]. In the DEM, the interaction between the contacting particles are treated as a dynamic process and the stress and deformation of the whole particles assembly are obtained from tracing the movement and interactions of each individual particle. The contact that connects any two particles can be represented through springs, friction resistance and damp absorber, as shown in Fig. 1 [22,23]. The dynamic process is completed through the integration of particles’ accelerations and velocities by using a central-difference scheme with an explicit time-step algorithm. The solution of a DEM model is carried out by applying Newton’s second law and the force–displacement law in turns to form one cycle of computation. Newton’s second law is used to calculate the particle’s acceleration resulting from the contact forces and external forces, while the force–displacement law is used to update the contact forces resulting from the overlap of the two contacting elements. These two laws are applied repeatedly to form the whole calculation cycle of DEM. Therefore, DEM is particular suitable to simulate the dynamic behavior of a particle system, in which the movement of every particle is recorded and analyzed over each time step. The discrete particles can also be densely packed and bonded together by adding special bonds at the contact points corresponding to special constitutive equations. This allows the DEM to simulate solid materials such as rock and concrete. The motion of each particle is determined by the resultant force and moment vectors acting upon it, and can be described in terms of translational and rotational motion [13,21]:

Dv i Dt Dwi Rotation Mi  bg wi ¼ I Dt

Translation F i  bg v i ¼ m

ð1Þ ð2Þ

where i (= 1, 2, 3) denotes, respectively the three co-ordinate directions, i.e., the x-, y-, z- directions; Fi is the out of balance force component of the particle; vi is the translational velocity; m is the mass of the particle; Mi is the out of balance moment due to contacts; wi is the rotational velocity; I is the rotational inertia of the particle; bg is the global damping coefficient and Dt is time step. For the 2D model used in the paper, the degrees of freedom of a particle element include the two motions in the x and y directions, respectively and the rotation about the z axis. For the force–displacement law, the resulting contact force vector at the contact between two particles, F, can be resolved into normal and shear components, Fn and Fs:

~ F ¼ F n~ n þ F s~ s

ð3Þ

The magnitude of the normal contact force is calculated by n

F ¼ K n Un

where Kn is normal stiffness at the contact and Un is the overlap of the two particles. The shear contact force is calculated by adding increment over a time step of Dt to the previous value, considering slippage is allowed to occur in the DEM:

F s ¼ F s þ DF s

ð5Þ

s

in which DF is calculated from the relative shear displacement DUs:

DF s ¼ K s DU s

ð6Þ

The contact is checked for slip conditions by calculating the maximum allowable shear contact force,

F s > ljF n j

ð7Þ

where l is the friction coefficient at the contact. 3. DEM modeling of fiber, matrix and interface 3.1. DEM modeling of fiber Solid materials are usually modeled by DEM through adding a bond at the contact of two contacting particles. Bonds in DEM can be envisioned as a kind of glue joining the two contacting particles. There are two intrinsic bonds, contact bond and parallel bond, in PFC2D which is a popular commercial code of DEM and will be used as the simulation platform of this research. A parallel bond can be regarded as a set of elastic springs with constant normal and shear stiffness, uniformly distributed over either a circular or rectangular cross-section lying on the contact plane and centered at the contact point, as shown in Fig. 2 [23]. Parallel bond can transmit both forces and moments. The contact bond can be viewed as a vanishingly small size of glue with a pair of elastic springs (or a point of glue) with constant normal and shear stiffness acting at the contact point. Therefore, contact bond can only transmit a force. The constitutive behavior of these two types of bonds is described in Fig. 3 [23]. It can be found in Fig. 3 that the contact bond and parallel bond have linear elastic behavior and fracture of the bonds is brittle, and thus they are suitable to model brittle materials. In this paper, parallel bond is chosen to represent the fiber material. In order to develop a DEM model for a single fiber, a virtual tension test was carried out to calibrate the contact properties between particles to benchmark the elastic properties obtained from laboratory tests. The stress–strain curve of the calibrated DEM model is shown in Fig. 4. All numerical simulations/tests in this paper were performed using PFC2D by running the programs written in FISH language, which is embedded in PFC2D enabling users to define new variables and functions. Random particles were packed to form a

ð4Þ

(a)

xi

(b)

A

ni

xi

C

xi B Particle is discal

Particle is spherical s

2R

n

Fi

Fi

n

Fi

s

Fi t

L L Fig. 1. Particles and contacts in DEM; (a) two particles in contact and (b) the physical elements of the contact.

Fig. 2. Parallel bond in DEM.

255

D. Yang et al. / Computational Materials Science 49 (2010) 253–259

Bond breaks

Fs

Fn Fcn

20

Bond breaks slip model

Fcs

v=0.02m/s v=0.05m/s v=0.1m/s v=0.2m/s v=1.0m/s

F smax

15 U

Us

(a) normal behavior

Stress (MPa)

slip model

n

(b) shear behavior

Fig. 3. Constitutive behavior of the bond at contact. (‘U’ is the overlap between the two particles, ‘F’ is the resulting contact force, ‘Fc’ is the bond strength, and superscripts ‘n’ and ‘s’ represent normal and shear direction, respectively.)

10

5

0 0.0

0.2

0.4

0.6

Strain (%) 150

σ

σ

Fig. 5. Convergence of loading velocities for quasi-static simulation.

Stress (MPa)

120 90 60 30 0 0.00

0.05

0.10

0.15

0.20

0.25

Strain (%) Fig. 4. Stress–strain curve in DEM simulation of fiber tensile test.

densely packed assembly, which was then subjected to a constant displacement rate at each ends. The fiber had a radius of 14.2 lm and a length of 1 mm, and the average particle radius was set to 1.25 lm. Since the length to average particle radius, L/R, is equal to 800, it is sufficiently large to ensure convergence of the simulation [16,24]. The details of the construction process of particle assemble have been described else where [16,19]. The numerical results are listed in Table 1. Since displacement loading was used in this simulation and the matrix tensile test and microbond tests shown in the following sections, the influence of loading velocities on the solutions was studied by comparing the stress–strain curves obtained for different loading velocities, as shown in Fig. 5. It is found that for all loading velocities below 0.1 m/s, the simulations lead to almost identical results. Therefore, v = 0.1 m/s will be chosen as the loading velocity in all the subsequent simulations to achieve the quasistatic conditions.

plastic deformation. Also, as suggested in [25,26] that the interface between fiber and matrix is adhesive, there exist residual interfacial traction forces, even when the fiber and matrix are detached but before they are entirely separated. Therefore alternative contact models have been proposed by DEM users to account for the complex interfacial behavior by considering more complicated constitutive laws. One of these is the contact softening model that is bilinear elastic and is based on the previous contact bond model [23,27]. The contact softening model is similar to the cohesive zone model (CZM) in the continuum mechanics [28,29]. The only difference between these two models is the unloading and reloading curves in the yielding part, as shown in Fig. 6. Contact softening model describes the behavior of contact bonds in elastic, and represents plastic deformation by linearly softening the bond after the contact force reaches the bond strength. In both tensile and shear situation, the bond strength decreased to zero when the plastic displacement reaches the maximum plastic displacement Upmax that is related to the fracture energy release rate G. The interfacial crack may behavior as mode I, mode II or mix mode according to the stress field in the crack tip. In order to automatically capture the three fracture models, the maximum plastic displacement Upmax was kept constant, while the bond normal and shear strengths were defined individually. Hence, in a 2D system, the fracture energy release rate for modes I and II are, respectively:

GI ¼

The epoxy matrix in composite system usually exhibits elastic– plastic behavior. The elastic bond models can rarely account for the Table 1 Material properties of the DEM model. Fiber modulus (glass) Fiber Poisson’s ratio Fiber radius Matrix modulus (epoxy) Matrix Poisson’s ratio Matrix yield stress Interface modulus Interfacial bond strength (tensile) Interfacial bond strength (shear) Interface fracture energy release rate (modes I or II)

60 GPa 0.23 14.2 lm 2.8 GPa 0.35 62 MPa 3.3 GPa 50 MPa 50 MPa 100 J/m2

ð8Þ

1  rs max  U p max 2

ð9Þ

and

GII ¼ 3.2. DEM modeling of matrix and interface

1  rn max  U p max 2

In 2D DEM, the contact stresses at the bond are taken as average stresses between the two contacting particles, as calculated below [23]:

Fn

Fs

Fcn

Fcs

U np

U npmax

(a) normal behavior

Un

U sp

U spmax

(b) shear behavior

Fig. 6. Constitutive behavior of contact softening model.

Us

256

D. Yang et al. / Computational Materials Science 49 (2010) 253–259

rk max ¼

F kc Fk ¼ c A 2Rt

ðk ¼ n; sÞ

1 R ¼ ðR1 þ R2 Þ 2

ð10Þ ð11Þ

where R1 and R2 are the radii of the two contacting particles, respectively; t is the dimension in the out-of-plane direction. In pure tensile or shear condition, the bond starts to yield when the contacting force exceeds the strength and breaks when the fracture energy release rate is reached. In the situation of mix mode fracture, the bond strength, Fmax, is calculated from the two strength parameters (i.e., F nc and F sc , in force unit) as a function of the current orientation of the contact force. It is assumed that the contact strength varies as linear function of the angle, a:

p

p

ð12Þ

where a is the angle between the directions of the contact force and the line segment connecting the centers of the two particles. The yielding of the bond in tension is determined by comparing the resultant contact force, i.e.,

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ¼ ðF n Þ2 þ ðF s Þ2

ð13Þ

with the contact strength. The bond yields if the contact force is larger than the contact strength:

F > F max

ð14Þ

The maximum plastic displacement Upmax is calculated from the displacements in the normal and shear directions. In each computation step, the increments of normal and shear plastic increments are subjected to a consistence condition of:

 2  2 ðDU p Þ2 ¼ DU np þ DU sp

ð15Þ

DU np DU sp

ð16Þ

¼ tan a

When the plastic displacement reaches Upmax,

X

jDU p j ! U p max

ð17Þ

The fracture energy release rate can be calculated as:

1 1  rn  U np þ  rs  U sp 2 2 X X 1 1 ¼  rn  jDU np j þ  rs  jDU sp j 2 2

G¼

ð18Þ

rn and rs are the normal and shear stress of the bond when yield occurs, i.e., Eqs. (13) and (14) are satisfied. For the mixed mode, the fracture energy release rate is somewhere between rates of the two pure fracture modes. Both the matrix and fiber/matrix interface are simulated by the contact softening model in this paper. Since this model will be used to stimulate fracture behavior, a convergence test to assess the effect of particle size has to be conducted first. A Double Cantilever Beam (DCB) virtual test [30,31], which has been widely applied to investigate interface strength, was conducted by 2D DEM model with cubic arrangement of particles to analyze tension failure of the contact softening model. The model consists of two layers made of the same material as shown in Fig. 7, each of which is further divided into two plies. The contacts between the particles within each layer were modeled by parallel bond, while the interface between the two layers was represented by contact softening model. The length of the plies is 40 mm and the thickness is 2 mm. An initial interfacial crack of 20 mm long from the left end of the

R=0.50mm R=0.25mm R=0.125mm R=0.0625mm

25

Loading Force (kN)

  2a 2a s  F nc þ F max ¼ 1   Fc

Fig. 7. DEM model of DCB test.

20 15 10 5 0 0

3

6

9

12

15

Loading Displacement (mm) Fig. 8. Convergence of particle sizes in the DEM simulations of DCB test.

model was generated by removing the interfacial bonds between the two layers. The right end of the specimen was fixed. The loading was applied by assigning a constant separation velocity in the vertical direction to the particles at left ends of the two layers, which results in a variable opening end displacement. Different particle sizes were used in this simulation to examine the particle size effect on the interfacial fracture and results are illustrated in Fig. 8. It is found that the model shows good convergence, though the results for R = 0.5 mm are a little scattered, which is due to the fact that only one ply of particle was used in each layer. The good convergence of the contact softening model indicates that the influence of particle size on the model is not significant and the model can be adopted to simulate the fracture of matrix and fiber/matrix interface. In Fig. 9, the contact force at the interface of the two layers against the opening displacements at the free ends is presented, where g = 0 represents the initial crack tip, as shown in the figure, the contact force gradually approaches zero at the tip (g = 0) when the end displacement increases, indicating that crack is propagating towards the interior zone of the beam. The contact force changes from tension to compression at a distance away from the crack tip. This analysis agrees very well the analysis of stress distribution of CZM by FEM and fracture mechanics with J-integration method [32]. The DEM model of epoxy matrix was calibrated through a tension test, in which the contacts of the particles were represented by the contact softening model, as shown in Fig. 10. In order to account for both the plastic deformation and final failure, bond failures are represented by colors1 of black and red assigned, respectively, when the bond starts to yield and when it finally breaks. It can be found that the material at both sides of the major crack has shown large amount of plastic deformation (concentration

1 For interpretation of color in Figs. 10, 14 and 16, the reader is referred to the web version of this article.

6

3 Fig. 12. DEM model of microbond test.

0

-3 0

1

2

3

4

5

6

Distance in front of crack tip η (mm) Fig. 9. Contact force distribution in front of the crack tip at different loading displacements, d (particle radius R = 0.25 mm).

σ

σ (a) initial model

(b) after failure Fig. 10. DEM model of epoxy tensile test. (Black dot indicates the bond has started to yield, and red dot indicates the final failure.)

60

Stress (MPa)

257

δ =4mm δ =6mm δ =8mm δ =10mm δ =11mm δ =12mm

9

n

Normal contact force Fc (kN)

D. Yang et al. / Computational Materials Science 49 (2010) 253–259

40

20

left end of the fiber and the vertical movement along the fiber center-line. The particles are between a minimum radius of Rmin = 1 lm and a maximum radius of Rmax = 1.5Rmin according to the uniform distribution law with a total particle number of 7618. The calibrated properties from Section 3 were assigned to the particles of different constituents and the contacts. The contacts between fiber particles (in green color)* and matrix particles (in blue color) were treated as the fiber/matrix interface. The vise blade was represented by a vertical rigid line with a much larger stiffness than those of the particles to avoid any deformation. And the blade was made to move horizontally with a constant velocity v = 0.1 m/s against the microbond specimen. A magnified view of the area near the vise blade is shown in Fig. 13, in which the contact network of the particles is plotted by connecting the centers of contacting particles with black lines. This network of contacts is similar to the mesh scheme in the FEM, but the contacts are substantially discrete and breakable. They can be assigned with different constitutive properties in order to simulate the fracture of multiphase materials. In this paper, two vise angles, u = 7° and u = 14°, were used in the DEM simulations, as shown in Fig. 14a and b. Complex pictures of the failure process during the microbond tests, such as dynamical damage initiation and propagation of matrix cracking as well as interfacial breaking, were all observed in the DEM modeling. DEM simulation results then were compared with experimental investigations which were shown in Fig. 15. An inclined vise with an angle of 30° was also used in the DEM simulation to exam the effect of the vise angle on the material damage, as shown in Fig. 16. The cutting forces related to the two vise angles were also compared and shown in Fig. 17. 5. Results and discussions

0 0

1

2

3

Strain (%)

In the cases of vertical vise, as shown in Fig. 14, the matrix in front of the vise, where the vise cut into, shows plastic deform (indicated by the concentration of black dots). The material closely surrounding the vise tip has been broken due to the high level of stress concentration. The red dots represent complete separation

Fig. 11. Stress–strain curve of epoxy matrix in tension.

of black dots). The stress–strain curve, as shown in Fig. 11, also demonstrates that the model has captured the elastic–plastic behavior of epoxy matrix accurately. The calibrated matrix properties are listed in Table 1. 4. DEM simulation of the microbond test In Fig. 12, a symmetric model of the microbond specimen was constructed by DEM, in which the droplet of epoxy matrix was modeled with an elliptical shaped particle assembly and the fiber center-line was taken as the symmetric axis. Configuration of the microbond test was taken from [9], i.e., the fiber radius was 14.17 lm, and the height and embedded length of the droplet were, respectively, 143.3 lm and 366.7 lm. The boundary conditions were applied by restricting the horizontal movement in the

Fig. 13. Network of the contacts in the DEM model of microbond test.

258

D. Yang et al. / Computational Materials Science 49 (2010) 253–259

(a) =7º

(b) =14º

(i)

(i)

(ii)

(ii)

(iii)

(iii)

Fig. 14. DEM simulation of microbond test at different vise angle u: (a) u = 7° and (b) u = 14°. (Black dot indicates the bond has started to yield, and red dot indicates the final failure.)

Matrix remaining on the fiber

Pulled-out droplet

Contact region with a knife edge Fig. 15. SEM experimental observation of microbond test in [10] and [33].

vertical vise, ϕ =7 vertical vise, ϕ =14 inclined vise, ϕ=7 inclined vise, ϕ=14

Cutting force (kN)

15

10

5

0 0

5

10

15

20

25

30

35

40

45

50

Cutting displacement (um) Fig. 17. Cutting force acting on the vise.

Fig. 16. DEM simulations of microbond test under inclined vise.

of the particles and breaking of bonds. The interface begins to yield before the matrix cracks propagate to it. The interface at the left end of the droplet has also yielded but debonding has not been observed throughout the microbond test. These observations agreed

D. Yang et al. / Computational Materials Science 49 (2010) 253–259

very well with the experimental investigations in Fig. 15 [10,33], shown also identical shape of the damaged matrix, the pattern of matrix fracturing and the appearance of plastic behavior on the pulled-out droplet. When the matrix cracking reached the interface, the interface debonding started and propagated along the interface, where plastic deformation within the matrix adjacent to the interface was observed, as shown in Fig. 14 (iii) and Fig. 15. When a larger vise angle (u = 14°) was used, as shown in Fig. 14b, more plastic deformation can be found in front of the vise, and matrix cracking tended to propagate towards the surface of droplet leading complete material detachment. This is due to the reduced contact area of matrix with the vise, and then less interface yielding to release the strain energy resulting from the shove of the vise. The larger the vise angle is, the shorter the length of interface debonding will be (Fig. 14 (iii)). When the vise angle is smaller, i.e., u = 7°, there is no significant difference in terms of material damage in the tests between using vertical and inclined vises, as shown in Fig. 16a. However, in the test of u = 14° (Fig. 16b), the matrix cracking did not propagate towards the surface due to the suppression of the inclined vise, and a larger plastic zone and more matrix cracking were induced underneath the vise tip. This also suggests that the material damage is sensitive to the vise angle when inclined vise is used. In other words, the geometry of the vise tip and the geometry of the droplet are both important in a microbond test, which to the authors’ best knowledge, has not well studied and documented in the literature [9,34]. The cutting forces (horizontal) acting on the vise were obtained from the DEM microbond tests and compared in Fig. 17. In the case of vertical vise, similar cutting forces were found in simulations under different vise angles. While in the DEM simulations of inclined vise, the cutting force was sensitive to the vise angle, and it was found that the larger the vise angle was used, the larger the maximum cutting force was induced. Furthermore, cutting forces in the simulations with inclined vise were found to be larger than those with vertical vises. Since the simulations were conducted in 2D, the calculation of interfacial shear strength was not included. 3D models have to be developed to better describe the post-failure friction behavior as well as include interfacial shear strength.

6. Conclusion A discrete element method (DEM) has been used to simulate the material failure process in the microbond test of fiber reinforced composites. The plastic deformation and the cracking of the matrix were simulated by introducing the contact softening model. The initiation and propagation of interfacial debonding were also captured naturally by the DEM simulation. Vertical and inclined vises with two different vise angles were used in the simulations. It was found that both the vise geometry and the vise angle had effects on the material damage process. It was also observed that the inclined vise was sensitive to the vise angle with a large variation of the cutting forces. The numerical results have been compared and validated by the existing experimental observations, which shows that the DEM model developed in the paper can be further modified to analyze material and structural damages in more complicated fiber reinforce composite systems.

259

Acknowledgements D.M. Yang would like to acknowledge the supports from the ORSAS scholarship and the University of Leeds studentship for his PhD study. Y.Q. Tan would like to thank the NSFC of China for the financial support of this collaboration research under the Grant No. 50875224. References [1] D. Hull, T.W. Clyne, An Introduction to Composite Materials, Cambridge University Press, 1996. [2] D. Trias, J. Costa, A. Turon, J.E. Hurtado, Acta Materialia 54 (2006) 3471–3484. [3] L.E. Asp, L.A. Berglund, R. Talreja, Composites Science and Technology 56 (1996) 657–665. [4] G. Lubineau, P. Ladeveze, Computational Materials Science 43 (2008) 137–145. [5] L. Mishnaevsky Jr., P. BrØdsted, Computational Materials Science 44 (2009) 1351–1359. [6] C.T. Chou, U. Gaur, B. Miller, Composites Science and Technology 51 (1994) 111–116. [7] X. Bi, Z. Li, P.H. Geubelle, J. Lambros, Mechanics of Materials 34 (2002) 433– 446. [8] N.-S. Choi, J.-E. Park, Composites Science and Technology 69 (2009) 1615– 1622. [9] J.T. Ash, W.M. Cross, D. Svalstad, J.J. Kellar, L. Kjerengtroen, Composites Science and Technology 63 (2003) 641–651. [10] M. Nishikawa, T. Okabe, K. Hemmi, N. Takeda, International Journal of Solids and Structures 45 (2008) 4098–4113. [11] R.J. Day, J.V.C. Rodrigez, Composites Science and Technology 58 (1998) 907– 914. [12] Itasca Consulting Group Inc., Minneapolis, Minnesota, 2004. [13] P.A. Cundall, O.D.L. Strack, Géotechnique 29 (1979) 47–65. [14] S.H. Frédéric, V. Donzé, L. Daudeville, Computers & Structures 82 (2004) 2509– 2524. [15] W.J. Shiu, F.V. Donze, L. Daudeville, Computers and Concrete 5 (2008) 329– 342. [16] D.O. Potyondy, P.A. Cundall, International Journal of Rock Mechanics and Mining Sciences 41 (2004) 1329–1364. [17] T.S. Wanne, R.P. Young, International Journal of Rock Mechanics and Mining Sciences 45 (2008) 789–799. [18] N. Cho, C.D. Martin, D.C. Sego, International Journal of Rock Mechanics and Mining Sciences 44 (2007) 997–1010. [19] Y. Tan, D. Yang, Y. Sheng, International Journal of Machine Tools and Manufacture 48 (2008) 975–982. [20] Y.Q. Tan, D.M. Yang, Y. Sheng, Journal of the European Ceramic Society 29 (2009) 1029–1037. [21] Y. Sheng, C.J. Lawrence, B. Briscoe, C. Thornton, Engineering Computations: International Journal for Computer-Aided Engineering 21 (2004) 304–317. [22] F.A. Tavarez, M.E. Plesha, International Journal for Numerical Methods in Engineering 70 (2007) 379–404. [23] Itasca Consulting Group Inc. (PFC2D user manual), Minneapolis, Minnesota, 2004. [24] B. Yang, Y. Jiao, S. Lei, Engineering Computations: International Journal for Computer-Aided Engineering 23 (2006) 607–631. [25] R. Borg, L. Nilsson, K. Simonsson, Composites Science and Technology 62 (2002) 1299–1314. [26] S. Li, M.D. Thouless, A.M. Waas, J.A. Schroeder, P.D. Zavattieri, Composites Science and Technology 65 (2005) 537–549. [27] H. Kim, M.P. Wagoner, W.G. Buttlar, Journal of Materials in Civil Engineering 20 (2008) 552–563. [28] D. Xie, A.M. Waas, Engineering Fracture Mechanics 73 (2006) 1783–1796. [29] M. Nishikawa, T. Okabe, N. Takeda, International Journal of Solids and Structures 44 (2007) 3101–3113. [30] M. Meo, E. Thieulot, Composite Structures 71 (2005) 429–434. [31] S. Li, M.D. Thouless, A.M. Waas, J.A. Schroeder, P.D. Zavattieri, Composites Science and Technology 65 (2005) 281–293. [32] Z.Y. Ouyang, G.Q. Li, Journal of Applied Mechanics-Transactions of the ASME, 76 (2009). [33] P. Zinck, H.D. Wagner, L. Salmon, J.F. Gerard, Polymer 42 (2001) 5401–5413. [34] B. Miller, U. Gaur, D.E. Hirt, Composites Science and Technology 42 (1991) 207–219.