New lattice models for dynamic fracture problems of anisotropic materials

Composites Part B 172 (2019) 760–768

Contents lists available at ScienceDirect

Composites Part B journal homepage: www.elsevier.com/locate/composites

New lattice models for dynamic fracture problems of anisotropic materials M. Braun b,c ,∗, M.P. Ariza a a

Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Sevilla 41092, Spain Departamento de Construcciones, Facultad de Ingeniería, Universidad Nacional de La Plata, Avda. 1 esq. 47, La Plata B1900TAG, Argentina c Consejo Nacional de Investigación Científicas y Técnicas (CONICET), CCT La Plata, Calle 8 No 1467, La Plata B1904CMC, Argentina b

ARTICLE

INFO

Keywords: Dynamic fracture Lattice model Composites Anisotropic materials Computational modeling Crack propagation

ABSTRACT In this paper we present a two-dimensional lattice model to simulate dynamic fracture problems in anisotropic materials. The model is based on the equivalence of strain energy stored in a representative unit cell to its counterpart continuum structure in the case of in-plane elasticity. A major advantage of this model is that the underlying lattice structure inherently incorporates the crack path preference of anisotropic materials. A numerical example is included to show the capability of these methods for modeling dynamic fracture problems in anisotropic materials. The results obtained are in agreement with the numerical results obtained with a finite element method. Discussions and future work are presented according to this research.

1. Introduction The use of fiber-reinforced composite materials has been proven to play a key role in various fields of engineering such as aerospace, automotive, marine and military engineering. The excellent strengthto-weight ratio of composite materials can be exploited in a variety of engineering applications for higher speed, longer range and greater payload capacity [1]. In recent decades, numerous numerical models have been developed to formulate and analyze progressive damage and failure in composite laminates. Continuum mechanics-based models include the finite element method (FEM) [2–5], the extended finite element method [6– 8], and meshless methods [9,10]. Although these models improve the accuracy and completeness of progressive damage analysis of composites, they still have limitations in dealing with discontinuities, such as dynamic crack propagation, multi-fragmentation and branching. There are other numerical methods based on discrete formulations, such as the discrete element method [11], peridynamics [12], and lattice models [13]. These models do not require an additional crack path determination rule and crack propagation is a natural outcome of the breakage of bonds. However, the development and application of discrete methods have mainly focused on isotropic materials. In this regard, Nayfeh and Hefzy [14] presented the first lattice model valid for an anisotropic material, but it had many limitations on the Young’s modulus selection. More recently, Xu et al. [15] used a peridynamic model to simulate the damage patterns in laminated composites subjected to bi-axial

loading and low-velocity impact. Their results in [16,17] show that homogenized peridynamic models are able to capture localization, quite well such as splitting fracture. In 2016, Chen et al. [18] developed a non-local lattice particle model for static fracture simulation of anisotropic materials. However, the fracture criteria of these models do not depend on the material strength properties. Therefore, it is not possible to apply classical composite failure criteria in these models. The main objective of this study was to extend the lattice particle model presented by Wang et al. [19] for isotropic materials to general anisotropic elastic solids and to explore the fracture behaviors of these types of materials. We considered two different lattice structures: square and triangular. The arbitrary material orientation of a general anisotropic material was taken into account by rotating the topological lattice structure. Additionally, we implemented the methodology proposed by Monette and Anderson [20] to calculate the stress tensor for lattice models (LM). Therefore, the proposed methods are the first discrete models capable of applying different composite failure criteria, such as the Tsai– Wu criterion [21], Tsai–Hill criterion [22] and Hashin criterion [23, 24]. In order to show the capability of the proposed LM, we studied the case of a rectangular composite plate submitted to a dynamic tensile load. Results are compared with a finite element analysis (FEA) implemented in Abaqus/Explicit [25].

∗ Corresponding author at: Departamento de Construcciones, Facultad de Ingeniería, Universidad Nacional de La Plata, Avda. 1 esq. 47, La Plata B1900TAG, Argentina. E-mail address: [email protected] (M. Braun).

https://doi.org/10.1016/j.compositesb.2019.05.082 Received 28 January 2019; Received in revised form 6 April 2019; Accepted 5 May 2019 Available online 7 May 2019 1359-8368/© 2019 Elsevier Ltd. All rights reserved.

Composites Part B 172 (2019) 760–768

M. Braun and M.P. Ariza

Table 1 Bond unit vectors n𝑏 in case of square particle interaction.

2. Description of the model 2.1. General considerations

𝑏

𝑏

𝐅 ⋅𝐮

where 𝑁 is the total number of bonds, 𝐅𝑏 is the axial force vector and 𝐮𝑏 is the displacement. The superscript 𝑏 is the 𝑏th bond. The force 𝐅𝑏 in the spring 𝑘𝑏 is related to the displacement 𝐮𝑏 by

8

𝑎2 ∑ 𝑘𝑏 𝐧𝑏𝑖 𝐧𝑏𝑗 𝐧𝑏𝑘 𝐧𝑏𝑚 4𝑎2 𝑏=1,3,5,7 √ ( 2𝑎)2 ∑ = 𝑘𝑏 𝐧𝑏𝑖 𝐧𝑏𝑗 𝐧𝑏𝑘 𝐧𝑏𝑚 4𝑎2 𝑏=2,4,6,8

I = 𝐶𝑖𝑗𝑘𝑚

(12)

II 𝐶𝑖𝑗𝑘𝑚

(13)

• Longitudinal direction: 𝐾1 = 𝑘1 = 𝑘5 • Transversal direction: 𝐾2 = 𝑘3 = 𝑘7 • Diagonal direction: 𝐾3 = 𝑘2 = 𝑘4 = 𝑘6 = 𝑘8

(4)

where 𝑎 is the initial bond distance. For the continuum model, the strain energy is:

substituting the Eq. (6) and the vectors presented in Table 1 into Eq. (11), we obtain:

𝑉 𝐶 𝜀 𝜀 (5) 2 𝑖𝑗𝑘𝑚 𝑖𝑗 𝑘𝑚 where 𝑉 is the area of the unit cell, 𝐶𝑖𝑗𝑘𝑚 is the fourth-order stiffness tensor, and 𝜀 is the linear strain tensor assumed to be uniform within the unit cell. Applying Eqs. (4) and (5) to Eq. (1), we obtain the stiffness tensor: 𝑎2 ∑ 𝑏 𝑏 𝑏 𝑏 𝑏 𝑘 𝐧𝑖 𝐧𝑗 𝐧𝑘 𝐧𝑚 (6) 𝐶𝑖𝑗𝑘𝑚 = 4𝑉 𝑏

𝑈𝑐 =

𝐸1 (1 − 𝜈12 𝜈21 ) 𝐸2 = (𝐾2 + 𝐾3 )∕2 = (1 − 𝜈12 𝜈21 ) 𝜈21 𝐸2 = 𝐾3 ∕2 = (1 − 𝜈12 𝜈21 )

𝐶1111 = (𝐾1 + 𝐾3 )∕2 =

(14)

𝐶2222

(15)

𝐶1122

𝐶1212 = 𝐾3 ∕2 = 𝐺12

(16) (17)

The stiffness of the lattice is obtained from the Eqs. (14)–(16):

The constants of the stiffness tensor as a function of engineering format in two-dimensional anisotropic material, considering plane stress, are described as: 𝐸1 𝐶1111 = (7) (1 − 𝜈12 𝜈21 ) 𝐸2 (8) 𝐶2222 = (1 − 𝜈12 𝜈21 ) 𝜈21 𝐸2 𝐶1122 = 𝐶2211 = (9) (1 − 𝜈12 𝜈21 )

2𝐸1 𝜈21 𝐸2 − (1 − 𝜈12 𝜈21 ) (1 − 𝜈12 𝜈21 ) 2𝐸2 1 𝐾2 = ( − 𝜈21 ) (1 − 𝜈12 𝜈21 ) 2 2𝜈21 𝐸2 𝐾3 = (1 − 𝜈12 𝜈21 ) 𝐾1 =

(18) (19) (20)

Thus, the constants of the lattice model are defined as a function of the material properties. We can define the Poisson’s ratios yields as follows:

(10)

𝐶1212 = 𝐺12

6 7

The model can be reduced to three different stiffness if we consider the symmetry of the unit vectors. Therefore, we can define the following relationships:

(3)

Therefore, the strain energy stored in a unit cell is: 𝑎2 ∑ 𝑏 𝑏 𝑏 𝑏 𝑏 𝑘 𝐧𝑖 𝐧𝑗 𝐧𝑘 𝐧𝑚 𝜀𝑖𝑗 𝜀𝑘𝑚 𝑈𝑑 = 8 𝑏

(0, −1) √ √ (−1∕ 2, −1∕ 2) (0, −1) √ √ (1∕ 2, −1∕ 2)

(2)

𝑏

𝐅 𝑏 = 𝑘𝑏 𝐮 𝑏

Diagonal: n𝑏

5

where

The discrete strain energy is given by the following equation: 1 𝑈𝑑 = 2

Diagonal: 𝑏

(1, 0) √ √ (1∕ 2, 1∕ 2) (0, 1) √ √ (−1∕ 2, 1∕ 2)

4

(1)

𝑁 ∑

Axial: n𝑏

1 2 3

The elastic properties relationship between the continuum and the lattice model is derived from the strain energy equivalence of both models [19,26]. We define the equivalence of discrete strain energy, 𝑈𝑑 , with the associated strain energy of the continuum, 𝑈𝑐 , 𝑈𝑑 = 𝑈𝑐

Axial: 𝑏

𝐾3 𝐶1122 = 𝐶2222 𝐾2 + 𝐾3 𝐾3 𝐸 = 2 𝐸1 (𝐾2 + 𝐾3 )

where 𝐸1 and 𝐸2 are the longitudinal and transverse elastic Young’s moduli in the principal material axes, respectively. 𝜈12 is the longitudinal Poisson’s ratio and 𝜈21 is the transverse Poisson’s ratio.

𝜈21 =

(21)

𝜈12

(22)

2.2. Square lattice model (SLM)

while the shear moduli is defined by the Eqs. (17) and (19)

Fig. 1 shows the square lattice shape and the square unit cell. For a particle 𝑖, there are eight axial and diagonal neighboring particles to be accounted for. A unit square cell area is 𝑉 = 𝑎2 . The bond vectors 𝐧𝑏 for each bond 𝑏 of a typical unit cell are given in Table 1. Note that, the axial bonds are listed in the left columns, whereas the diagonal bonds are listed in the right ones. It should be noted that in this lattice shape each particle communicates with four neighbors (𝑏 = 1, 3, 5, 7) with a distance equal to 𝑎, while the√distance between neighbors in the diagonal (𝑏 = 2, 4, 6, 8) is equal to 2𝑎. If we consider the equivalence of strain energy in a unit cell and an effective continuum, we obtain the following effective local-type stiffness tensor [26]:

𝐺12 =

𝐶𝑖𝑗𝑘𝑚

= 𝐶I

𝑖𝑗𝑘𝑚

+ 𝐶 II

𝑖𝑗𝑘𝑚

𝜈21 𝐸2 (1 − 𝜈12 𝜈21 )

(23)

2.3. Triangular lattice model (TLM) The triangular lattice structure and the hexagonal unit cell are √ shown in Fig. 2. The unit cell area is 𝑉 = 3𝑎2 ∕2, where 𝑎 is the distance between particles. For each bond 𝑏, 𝑘𝑏 is the spring constant of half-length of the bond, and the unit cell bond vectors 𝐧𝑏 are given in Table 2. We considered only two different stiffness values due to the symmetry of the unit cell bond vectors. • Longitudinal direction: 𝐾1 = 𝑘1 = 𝑘4 • Transversal direction: 𝐾2 = 𝑘2 = 𝑘3 = 𝑘5 = 𝑘6

(11) 761

Composites Part B 172 (2019) 760–768

M. Braun and M.P. Ariza

Fig. 1. A square lattice with a square unit cell.

Fig. 2. A triangular lattice with a hexagonal unit cell. Table 2 Bond unit vectors n𝑏 in case of triangular particle interaction.

From Eqs. (24) and (25), we obtain the stiffness as a functions of the mechanical material properties:

𝑏

𝐧𝑏

𝑏

𝐧𝑏

1

(1, 0)

4

(−1, 0)

2 3

√ (1∕2, 3∕2) √ (−1∕2, 3∕2)

5 6

8𝐸2 𝐾2 = √ 3 3(1 − 𝜈12 𝜈21 ) √ 3𝐸1 8𝐸2 𝐾1 = − √ (1 − 𝜈12 𝜈21 ) 3 3(1 − 𝜈 𝜈 ) 12 21

√ (−1∕2, − 3∕2) √ (1∕2, − 3∕2)

𝜈21 𝐸2 1 3𝐾2 𝐶1122 = √ = 4 (1 − 𝜈12 𝜈21 ) 2 3 1 3𝐾2 𝐶1212 = √ = 𝐺12 2 3 4

(29)

Applying the other three equations, we obtain the following limitations: 𝐸 𝐺12 = 2 (30) 3 𝐶1122 1 𝜈21 = = (31) 𝐶2222 3 𝐸 𝜈12 = 2 (32) 3𝐸1

From Eq. (6), we can derive the following: 𝐾 𝐸1 1 𝐶1111 = √ (2𝐾1 + 2 ) = 4 (1 − 𝜈12 𝜈21 ) 2 3 9𝐾 𝐸 1 2 2 𝐶2222 = √ = (1 − 𝜈12 𝜈21 ) 2 3 4

(28)

(24) (25) (26)

Note that 𝐸2 ≪ 𝐸1 , because direction 1 corresponds to the fiber direction. Therefore, 𝜈12 has values between 0 and 0.5. Limitations in the choice of Poisson’s ratios are also present in other lattice models [19, 27–29].

(27)

762

Composites Part B 172 (2019) 760–768

M. Braun and M.P. Ariza

The equation for temporal advancing positions is 𝐮(𝑡 + 𝛥𝑡) = 2𝐮(𝑡) − 𝐮(𝑡 − 𝛥𝑡) + 𝛥𝑡2 𝐮̈ (𝑡) + 𝑂(𝛥𝑡4 ).

(33)

To ensure that numerical errors did not grow rapidly in time, we applied the Courant–Friedrichs–Lewy criterion [31], according to which the time increment depends on the discrete spacing 𝑎 and the wave speed 𝑐𝑝 , 𝛥𝑡 ≤

𝑎 . 𝑐𝑝

(34)

where 𝑐𝑝 is calculated as the maximum value between √ 𝐸2 ∕𝜌, while 𝜌 is the material density.

√ 𝐸1 ∕𝜌 and

3. Failure behavior 3.1. Stress computation The stress acting upon a cubic cell are usually defined as 𝜎𝑘𝑙 =

∑

𝐅𝑏

𝑏

𝐧𝑏𝑘𝑙

(35)

𝐴

where 𝐅𝑏 is the force on any surface 𝑏 of the cubic cell, while 𝐧𝑏𝑘𝑙 is a unit vector (𝑘𝑙 = 11, 22, 33) or parallel (𝑘𝑙 = 12, 23, 13) to the surface 𝑏 and 𝐴 is the surface area. This definition was adapted to two-dimensional triangular and square lattices by [20] as follows: (𝜎11 )𝑖 = (𝜎22 )𝑖 = (𝜎12 )𝑖 = (𝜎21 )𝑖 =

𝑚 ∑ 𝑏=1 𝑚 ∑ 𝑏=1 𝑚 ∑ 𝑏=1 𝑚 ∑ 𝑏=1

Fig. 3. Correspondence between material orientation and lattice rotation.

(𝐅𝑏 )1 (𝐅𝑏 )2 (𝐅𝑏 )1 (𝐅𝑏 )2

(𝐧𝑏𝑘𝑙 )1 𝑙0 (𝐧𝑏𝑘𝑙 )2 𝑙0 (𝐧𝑏𝑘𝑙 )2 𝑙0 (𝐧𝑏𝑘𝑙 )1 𝑙0

(36) (37) (38) (39)

∑ where 𝑚 𝑏=1 is a sum over 𝑚 neighbors of particle 𝑖, (𝐅𝑏 )𝑛 is the force component in the direction 𝑥𝑛 acting on node 𝑖 from node 𝑏, and (𝐧𝑏𝑘𝑙 )𝑛 is the component in the direction 𝑥𝑛 of √ the unit vector to the cell surface. Finally, 𝑙0 is a surface length, 𝑙0 = 3𝑎∕2 for the triangular lattice and 𝑙0 = 𝑎 for the square lattice.

Moreover, we were only able to effectively match parameter 𝐸1 and 𝐸2 for the unidirectional composite. The other material parameters (Poisson’s ratios, shear moduli) depend on these two values. However, the model presented by Hu et al. [16] has the same limitation. Nevertheless, the shear moduli and Poisson’s ratios have a secondary importance in the effect of the dynamic fracture and damage in unidirectional materials [16].

3.2. Failure criterion There are various discrete failure criteria for isotropic materials in the literature, such as: critical force [32], critical strain energy [33], maximum principal stress [34–36], and critical elongation [37]. Discrete failure criteria developed for anisotropic materials [16,18] are based on a critical elongation. This parameter is not based on a physical property, and requires a calibration procedure to define it. However, the discrete model developed in this study is capable of incorporating the classical failure criteria for composite materials [38] (e.g., Tsai–Hill, Hashin, Tsai–Wu). This was possible because we proposed a methodology to calculate the stress tensor for each bond (see Section 3.1). In this study, we implemented Hashin’s theory [23,24] to compare our model with the finite element method implemented in Abaqus/Explicit commercial software [25]. This criterion consider four different conditions:

2.4. Material orientation The definition of the lattice orientation and how to rotate the lattice according to the material orientation is schematically shown in Fig. 3. The change in the fiber orientation consists of rotating the initial lattice (triangular or square) in a certain angle with respect to the horizontal axis, generating the geometry of the plate according to this new orientation. The lattice orientation always coincides with the material coordinate orientation and can be represented by the underlying lattice orientation. As shown in Fig. 3, the underlying lattice is rotated according to the rotation of the material coordinate following the same procedure presented by Chen et al. [18].

• Fiber tension ((𝜎11 )𝑖 > 0): ( ) ) ( (𝜎11 )𝑖 2 (𝜎12 )𝑖 2 + 𝛼 ⩾1 𝑋𝑇 𝑆𝐿

2.5. Time integration We resorted to the Verlet algorithm [30] for the time integration, which provides a direct solution for second-order differential equations. This method does not involve velocities calculation, reducing computational costs, but it is based on positions 𝐮(𝑡), accelerations 𝐮̈ (𝑡), and positions from the previous step 𝐮(𝑡 − 𝛥𝑡).

• Fiber compression ((𝜎11 )𝑖 < 0): ( ) (𝜎11 )𝑖 2 ⩾1 𝑋𝐶 763

(40)

(41)

Composites Part B 172 (2019) 760–768

M. Braun and M.P. Ariza Table 3 Mechanical properties of the plate. Property Longitudinal Young’s modulus 𝐸1 Transversal Young’s modulus 𝐸2 Shear modulus 𝐺12 Longitudinal Poisson’s ratio 𝜈12 Transversal Poisson’s ratio 𝜈21 Density 𝜌 Longitudinal Fracture Energy 𝐺1 Transversal Fracture Energy 𝐺2 Fiber strengths 𝑋 𝐶 = 𝑋 𝑇 Matrix strengths 𝑌 𝐶 = 𝑌 𝑇 Shear strengths 𝑆 𝐿 = 𝑆 𝑇 Coefficient 𝛼

134 GPa 7 GPa 4.4 GPa 0.33 0.0825 1530 kg/m3 15.49 kJ/m2 0.168 kJ/m2 1270 MPa 42 MPa 93 MPa 0

Fig. 4. Geometry and boundary conditions of a rectangular plate with a center hole.

• Matrix tension ((𝜎22 )𝑖 > 0): ( ) ( ) (𝜎22 )𝑖 2 (𝜎12 )𝑖 2 + ⩾1 𝑇 𝐿 𝑌 𝑆 • Matrix compression ((𝜎22 )𝑖 < 0): [( ] ( ) )2 ( ) (𝜎22 )𝑖 (𝜎22 )𝑖 2 (𝜎12 )𝑖 2 𝑌𝐶 + − 1 + ⩾1 2𝑆 𝑇 2𝑆 𝑇 𝑌𝐶 𝑆𝐿

(42)

Fig. 5. Time evolution of SLM reaction load for different mesh sizes.

of the plate were 200 × 100 × 2.5 mm3 [16], the plate had a hole of 25 mm in diameter, as detailed in Fig. 4. The composite material adopted in this case corresponds to a M55J/M18 carbon/epoxy. All mechanical properties of the plate are shown on Table 3. A traction velocity 𝑣 of 2 m/s was applied on the top face, while the lower face was fixed and the lateral faces are free. The velocity was applied as a Heaviside function starting at 𝑡 = 0.

(43)

where 𝑋, 𝑌 are the lamina longitudinal and transverse strengths, while the indices 𝑇 and 𝐶 correspond to the tension and compression resistance, respectively. 𝑆 𝐿 and 𝑆 𝑇 denote the longitudinal and transverse shear strength. 𝛼 is a coefficient that determines the contribution of the shear stress to the fiber tensile initiation criterion, while (𝜎11 )𝑖 , (𝜎22 )𝑖 and (𝜎12 )𝑖 are the lamina longitudinal, transversal and shear stress of the particle 𝑖, respectively. The stress components are calculated from Eqs. (36)–(39). The formulation of the fracture criterion is as follows: when any condition of Hashin’s criterion is not verified, the interaction between the node and its neighbors is assumed to be deleted, and the crack propagates through the created void. Note that, the interactions are removed when the toughness of the bonds is zero. This formulation is based on the methodology implemented in [39] for isotropic materials.

4.2. Mesh size analysis In this section, we analyze the mesh size dependence of LM proposed in this research. We analyzed the case detailed in Section 4.1 considering a fiber angle inclination of 45◦ . We explored the effect of five mesh sizes: 5.0, 4.0, 3.5, 2.5 and 1.25 mm. Figs. 5 and 6 present the time evolution of the reaction force obtained with SLM and TLM for different mesh sizes. It can be observed that the proposed lattice models converge on a same solution in the elastic range for all mesh sizes, in both models. Moreover, the mesh size influenced the results after the peak load, when crack propagation began. Both models showed the same trend: when mesh size decreased, the force load curve tended to move to the left. In addition, the difference between the results decreased when the mesh size was reduced. Figs. 7 and 8 show the crack pattern obtained for different levels of discretization of SLM and TLM, respectively. The crack inclination angle was approximately 45◦ for TLM, in agreement with the fiber inclination angle for all mesh sizes. Therefore, the crack pattern obtained with TLM is defined by the matrix tension failure. The crack angle obtained with SLM was lower than the angle obtained with TLM. In this case, the crack propagation was defined by the fiber and the matrix tension failure. The crack pattern did not

4. Numerical validation In this section, we explore the capability of the proposed models to capture the main features of the crack propagation process, determined by the fracture load and crack patterns. We considered the case of a rectangular composite plate with a center hole under imposed traction velocity. We compared the lattice model results with an FEA to validate the proposed methods. 4.1. Problem description The fracture process of a unidirectional composite plate with a center hole under imposed traction velocity was modeled. The dimensions 764

Composites Part B 172 (2019) 760–768

M. Braun and M.P. Ariza

Table 4 CPU time in function of the element size obtained with SLM and TLM. Size [mm]

SLM [s]

TLM [s]

5.00 4.00 3.50 2.50 1.25

73.44 136.45 188.08 563.64 9620.88

74.15 119.42 176.77 416.29 5410.58

problem with intermediate mesh sizes, the computational cost is not an important factor to choose between the SLM or TLM. The relatively low computational costs prove the simplicity of the model. We showed that the results for 𝑎 = 2.50 mm had a reasonable precision. We also proved that the proposed models converged into one solution. Fig. 6. Time evolution of TLM reaction load for different mesh sizes.

4.3. Fiber orientation dependence have a strong dependence on the mesh size but rather changed with the discretization type. Table 4 presents the values of CPU times taken to solve the problem, for each size mesh obtained with SLM and TLM. As we expected, for the smaller mesh sizes CPU time was greater in the SLM than in the TLM. The reason for this difference is that each particle 𝑖 of the SLM had eight neighbors while the particles of the TLM have six. Note that the CPU time of the SLM and TLM were approximately the same for mesh sizes greater than 2.50 mm. Therefore, if we explore a

To explore the capability of the proposed models for capturing the fiber orientation dependence, we considered six inclinations angles: 0, 30, 45, 60 and 90 degrees. Based on the results presented in the last section, we justified the discretization level adopted in the simulation, that is, 𝑎 = 2.5 mm for both lattice structures. Results obtained with the SLM and TLM were compared with the FEM. The FEA was implemented using Abaqus/Explicit commercial software [25]. The plate was meshed using S4R elements (a 4-node

Fig. 7. Crack pattern obtained with SLM for different mesh sizes.

Fig. 8. Crack pattern obtained with TLM for different mesh sizes. 765

Composites Part B 172 (2019) 760–768

M. Braun and M.P. Ariza

Fig. 9. Time evolution of the reaction force for different 0◦ fiber orientation angle.

Fig. 11. Time evolution of the reaction force for different 45◦ fiber orientation angle.

Fig. 10. Time evolution of the reaction force for different 30◦ fiber orientation angle.

Fig. 12. Time evolution of the reaction force for different 60◦ fiber orientation angle.

doubly curved thin or thick shell, reduced integration, finite membrane strains). It was discretized with 2979 elements (the element sizes are 2.5 mm). Damage behavior was defined by Hashin’s theory, implemented in Abaqus/Explicit software [25]. Figs. 9–13 show the reaction force as a function of the time calculated for different fiber orientation angles. The square and triangular lattice models had a good agreement with the FEM in the elastic range. However, the TLM had a smaller difference with the FEM than SLM in all cases. Note that the lattice models are valid for small strains, so the difference with the FEM increases with time. After the peak loads, the difference between the lattice models and the FEM increased since the proposed models did not consider the fracture energies, while the FEA included a softening law as a function of those energies. This difference increased when the fracture process was controlled by the fiber tension, because the fiber fracture energy was greater than the matrix fracture energy.

Fig. 13. Time evolution of the reaction force for different 90◦ fiber orientation angle.

Table 5 shows the failure load values obtained with the FEM, SLM and TLM. Moreover, there were relative errors of the SLM and TLM compared to the FEM. The results obtained show the capability of the proposed methods for calculating the stress tensor. Furthermore, the LM was in agreement with the FEM, showing the maximum difference at 0 and 90 degrees for the SLM and TLM, respectively. The upper SLM error was 15.2%, while the upper TLM error was 10.8%. The lowest difference was obtained for a 45◦ fiber angle.

Figs. 14–18 present the crack pattern obtained with the three methods for each fiber angle analyzed. The crack trajectory obtained with the FEM is represented by the red elements. We can see that the crack inclination angle obtained with FEM did not change for the fiber angles lower than 60◦ , while the crack propagation angle of LM increased with the fiber angle inclination. The 90◦ fiber angle exhibited a complex crack pattern since the fracture was defined by a mixed failure between the matrix and the fiber tension.

Therefore, the proposed models were able to calculate the failure load using Hashin’s failure criterion. 766

Composites Part B 172 (2019) 760–768

M. Braun and M.P. Ariza Table 5 Failure load of a plate with a hole obtained for each fiber orientation angle. Fiber angle [degrees]

FEM [kN]

SLM [kN]

TLM [kN]

SLM [relative error]

TLM [relative error]

0 30 45 60 90

5.89 6.59 8.13 11.12 104.69

6.79 7.07 8.37 11.95 113.33

5.81 6.21 8.18 10.23 93.43

15.2% 7.3% 2.9% 7.5% 8.3%

1.3% 5.7% 0.6% 8.0% 10.8%

Fig. 16. Crack pattern obtained for 45◦ fiber orientation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 14. Crack pattern obtained for 0◦ fiber orientation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 17. Crack pattern obtained for 60◦ fiber orientation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 15. Crack pattern obtained for 30◦ fiber orientation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The crack patterns obtained based on the proposed model are in agree with experimental data for low strain rate responses [40,41]. These experimental results show that the crack angle is parallel to the fiber direction, for low strain rates. Results show that the triangular lattice was more mesh-dependent than the square lattice. However, this dependence could be an advantage for anisotropic materials where the crack trajectories have preferential ways. Results show the good agreement between the proposed methods with the FEM but do not make it possible to determine which results are best. Thus, our models need to be compared with experimental results for composite materials submitted to dynamic loads.

Fig. 18. Crack pattern obtained for 90◦ fiber orientation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

anisotropic materials. The proposed models were based on the equivalence of strain energy between the discrete and the continuum systems. Lattice particles interact only with the nearest neighbors. Interactions were defined by axial bonds that depended on the stiffness matrix. Moreover, we adopted the methodology proposed by Monette and Anderson [20] to calculate the stress tensor for each particle. The

5. Conclusions In this research, we proposed two-dimensional square and triangular lattice models for the simulations of dynamic fracture problems of 767

Composites Part B 172 (2019) 760–768

M. Braun and M.P. Ariza

proposed methodology is the first discrete method that is able to apply the classical composite failure criteria. In order to evaluate the applicability of the proposed lattice, we modeled a rectangular composite plate with a center hole submitted to dynamic tensile load. Results were compared to finite element solutions, implemented in commercial software. The models developed in this study were able to obtain the same solution quality of the FEM when applied to dynamic crack propagation problems in unidirectional materials. Future studies are required to further quantitatively validate the model against experimental results for composite materials. Furthermore, we can study the random nature of the material using an aleatory distribution of mechanical properties associated with each bond. Additional research is also required to extend the proposed model to a general three-dimensional case.

[14] Nayfeh AH, Hefzy MS. Continuum modeling of three-dimensional truss-like space structures. AIAA J 1978;16(8):779–87. [15] Xu J, Askari A, Weckner O, Silling S. Peridynamic analysis of impact damage in composite laminates. J Aerosp Eng 2008;21(3):187–94. [16] Hu W, Ha YD, Bobaru F. Peridynamic model for dynamic fracture in unidirectional fiber-reinforced composites. Comput Methods Appl Mech Engrg 2012;217:247–61. [17] Hu W, Ha YD, Bobaru F. Modeling dynamic fracture and damage in a fiberreinforced composite lamina with peridynamics. Int J Multiscale Comput Eng 2011;9(6). [18] Chen H, Jiao Y, Liu Y. A nonlocal lattice particle model for fracture simulation of anisotropic materials. Composites B 2016;90:141–51. [19] Wang G, Al-Ostaz A, Cheng A-D, Mantena P. Hybrid lattice particle modeling: theoretical considerations for a 2d elastic spring network for dynamic fracture simulations. Comput Mater Sci 2009;44(4):1126–34. [20] Monette L, Anderson M. Elastic and fracture properties of the two-dimensional triangular and square lattices. Modelling Simulation Mater Sci Eng 1994;2(1):53. [21] Tsai SW, Wu EM. A general theory of strength for anisotropic materials. J Compos Mater 1971;5(1):58–80. [22] Azzi V, Tsai S. Anisotropic strength of composites. Exp Mech 1965;5(9):283–8. [23] Hashin Z. Failure criteria for unidirectional fiber composites. J Appl Mech 1980;47(2):329–34. [24] Hashin Z, Rotem A. A fatigue failure criterion for fiber reinforced materials. J Compos Mater 1973;7(4):448–64. [25] Hibbitt K, Sorensen. Abaqus explicit v6.13 user’s manual, version 6.13. Richmond, USA: Abaqus Inc.; 2013. [26] Grah M, Alzebdeh K, Sheng P, Vaudin M, Bowman K, Ostoja-Starzewski M. Brittle intergranular failure in 2d microstructures: experiments and computer simulations. Acta Mater 1996;44(10):4003–18. [27] Martín T, Espanol P, Rubio MA, Zúniga I. Dynamic fracture in a discrete model of a brittle elastic solid. Phys Rev E 2000;61(6):6120. [28] Wang G, Al-Ostaz A, Cheng A-D, Mantena P. A macroscopic-level hybrid lattice particle modeling of mode-i crack propagation in inelastic materials with varying ductility. Int J Solids Struct 2009;46(22–23):4054–63. [29] Kosteski L, DAmbra RB, Iturrioz I. Crack propagation in elastic solids using the truss-like discrete element method. Int J Fracture 2012;174(2):139–61. [30] Allen M, Tildesley D. Computer simulation of liquids. New York: Oxford University Press; 1987. [31] Bathe K. Finite element procedures. New Jersey: Prentice Hall; 1996. [32] Jagota A, Bennison S. Element breaking rules in computational models for brittle fracture. Modelling Simulation Mater Sci Eng 1995;3(4):485. [33] Silling S, Weckner O, Askari E, Bobaru F. Crack nucleation in a peridynamic solid. Int J Fract 2010;162(1–2):219–27. [34] Braun M, Fernández-Sáez J. A new 2d discrete model applied to dynamic crack propagation in brittle materials. Int J Solids Struct 2014;51(21–22):3787–97. [35] Braun M, Fernández-Sáez J. A 2d discrete model with a bilinear softening constitutive law applied to dynamic crack propagation problems. Int J Fract 2016;197(1):81–97. [36] Braun M, González Albuixech VF. Analysis of the stress intensity factor dependence with the crack velocity using a lattice model. Fatigue Fract Eng Mater Struct 2019;42(5):1075–84. [37] Bolander JE, Sukumar N. Irregular lattice model for quasistatic crack propagation. Phys Rev B 2005;71(9):094106. [38] Barbero EJ. Introduction to composite materials design. CRC press; 2017. [39] Martin T, Espanol P, Rubio M. Mechanisms for dynamic crack branching in brittle elastic solids: strain field kinematics and reflected surface waves. Phys Rev E 2005;71(3):036202. [40] Haque A, Ali M. High strain rate responses and failure analysis in polymer matrix composites - An experimental and finite element study. J Compos Mater 2005;39(5):423–50. [41] Li X, Yan Y, Guo L, Xu C. Effect of strain rate on the mechanical properties of carbon/epoxy composites under quasi-static and dynamic loadings. Polym Test 2016;52:254–64.

Acknowledgments The authors are indebted to the AUIP (Asociación Universitaria Iberoamericana de Postgrado), Spain and the Consejería de Economía, Innovación, Ciencia y Empleo of Junta de Andalucía, Spain (P12-TEP850) for the financial support. The authors express sincere gratitude to Professor José Fernández-Sáez for helpful discussions. References [1] Vinson J. The behavior of sandwich structures of isotropic and composite materials. Routledge; 2018. [2] Chang F-K, Chang K-Y. A progressive damage model for laminated composites containing stress concentrations. J Compos Mater 1987;21(9):834–55. [3] Chang K-Y, Llu S, Chang F-K. Damage tolerance of laminated composites containing an open hole and subjected to tensile loadings. J Compos Mater 1991;25(3):274–301. [4] Hallett S, Green B, Jiang W, Wisnom M. An experimental and numerical investigation into the damage mechanisms in notched composites. Composites A 2009;40(5):613–24. [5] Moure M, Sanchez-Saez S, Barbero E, Barbero E. Analysis of damage localization in composite laminates using a discrete damage model. Composites B 2014;66:224–32. [6] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Internat J Numer Methods Engrg 1999;45(5):601–20. [7] Wang Y, Waisman H. Progressive delamination analysis of composite materials using XFEM and a discrete damage zone model. Comput Mech 2015;55(1):1–26. [8] Yazdani S, Rust WJ, Wriggers P. An XFEM approach for modelling delamination in composite laminates. Compos Struct 2016;135:353–64. [9] Nguyen VP, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79(3):763–813. [10] Liew KM, Zhao X, Ferreira AJ. A review of meshless methods for laminated and functionally graded plates and shells. Compos Struct 2011;93(8):2031–41. [11] Maheo L, Dau F, Andre D, Charles J-L, Iordanoff I. A promising way to model cracks in composite using discrete element method. Composites B 2015;71:193–202. [12] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 2000;48(1):175–209. [13] Van Mier JG, van Vliet MR, Wang TK. Fracture mechanisms in particle composites: statistical aspects in lattice type analysis. Mech Mater 2002;34(11):705–24.

768