Study of Dynamic Brittle Fracture of Composite Lamina Using a Bond-Based Peridynamic Lattice Model

We proposed a bond-based peridynamic lattice model for simulating dynamic brittle fracture of 2D composite lamina. Material orthogonal anisotropy was represented by rotating topological lattice structure instead of fiber directions. Analytical derivation and numerical implementation of the proposed model were given based on energy equivalence. Benchmark composite lamina tests are used to validate the capability of modeling dynamic fracture of the method. )e peridynamic lattice model is found to be robust and successful inmodeling dynamic brittle fracture of 2D composite lamina and can be extended to composite laminates by applying 3D lattice structure.


Introduction
Fiber-reinforced composite materials [1,2] have been widely used in many primary structural components in various fields, such as aerospace, naval architecture, and automobile, due to their outstanding material properties like high strength-to-weight ratio and high stiffness-to-density ratio.Composite materials are classical inhomogeneous, anisotropic, and brittle materials, and understandings of their failure mechanisms and failure modes are still challenges for current mainstream computational methods, like finite element method (FEM) [3,4].
FEM is a method based on continuum mechanics, which has several limitations in modeling discontinuities.Firstly, the stress/strain conception is not applied where discontinuities emerge, and damage onset and propagation cannot be defined without extra displacement treatments; secondly, even though discontinuities can be dealt with some special techniques, like redefining a body to exclude the crack and then applying conditions at the crack surfaces as boundary conditions, mesh refine technique is still needed to obtain robust and accurate results; thirdly, the FEM model suffers from the requirement of extra damage criterion of guiding how crack propagates.All above limitations are caused by the differential form governing equations of motion of classical continuum mechanics.
In order to solve these limitations, Silling [5]and Silling et al. [6,7] introduced a nonlocal theory which uses spatial integration instead of its derivatives, called peridynamics (PD).In PD models, internal forces are represented by integration of nonlocal interactions between material points within a family, which is defined by a nonlocal character length, and damage information is involved in the material constitutes.With advantages mentioned above, due to the nonlocal characteristic, PD models suffer from larger computation cost and lower computation efficiency than that of FEM.Besides, the inhomogeneity of composite materials leads to a larger number of nodes and a finer mesh.
Bond-based PD has already been successfully applied in the prediction failure of composites.Askari et al. [8] and Xu et al. [9] predicted the strength of composite with notches under biaxial loadings.Colavito et al. [10] predicted damage propagation of unidirectional reinforced laminated plate under impact loading and static indentation.Hu et al. [11] proposed a model for composite lamina fracture by bridging the microscale peridynamic bonds and the macroscale parameters of the lamina; Hu et al. [12] developed a model for simulating delamination by introducing an energy-based approach, which removes bonds between adjacent layers.
However, low computation efficiency hinders applications of the PD model, for example, Oterkus et al. [13] combined FEM with peridynamics in predicting failure of a large composite panel with a central slot; still, large amounts of computation resources are consumed; thus, it is desirable to seek methods of improving efficiency.Besides coupling with FEM, other practical methods have been proposed, for example, parallel implementation [14] by Killic et al., homogeneous method [15] by A. Askari et al., domain decomposition method applied by Killic et al., and harmonic function for representing fiber direction in a uniform way by Askari et al. [15] and Azdoud et al. [16] to some extent reduces computation cost; however, the computation cost is far beyond the ability of personal computers, and thus, a practical solution should be developed.As pointed by Gerstle [17] and Gerstle [18], for many problems in solid mechanics which do not deform too much during a simulation, it can be assumed that material points (particles) will not change their neighbors.en, computation savings can be obtained by predefining neighbor of particles in a fixed topology before entering time integration loop.
e predefined neighborhood avoids a decomposition of the body domain, and this can be done by using the lattice structure.
Lattice modeling technique, first proposed by Hrennikoff [19], has been successfully applied in simulation of elastic and dynamic fracture of brittle materials, such as ceramics, metals, and composites [20][21][22], and high computation efficiency was achieved.However, mathematical structure of lattice methods is of Eulerian; its applications are limited to static fractures and sometimes microstructure systems.To take advantage of high computation efficiency and direction sensitive of the lattice model and extend its application to dynamic and macrolevel problems, it is suitable to transplant a lattice model into a new dynamic model, for example, a hybrid lattice particle model by Wang [23] and the state-based peridynamic lattice model (SPLM) by Gerstle et al. [17,18], which has been applied to the strength and failure prediction of concrete structures.However, they did not explore the influence of lattice structure and extend the method to composite materials.
In this paper, a hybrid method of (bond-based) peridynamic model and lattice model (PDLM) was proposed to simulate fracture of composite lamina.In order to model the orthotropic property of composite lamina, bonds (springs) of hexagonal lattice are assigned to different properties: bonds parallel to fiber direction are defined as fiber bonds, and all other directions are matrix bonds, and the assignment can be done using homogeneous technique from [14] by introducing a spherical function.To simulate fracture process of composite lamina under uniaxial loading, external loads are imposed on two edges of composite plate horizontally.e special boundary treatment for body force and constraints are not needed anymore, and a simple boundary displacement assumption was made.For the elastic case, bonds are assumed to be unbroken, and displacement contour is compared to that of analytical; for the brittle fracture case, bond breaking is controlled by a critical stretch, and when a bond breaks, force in the bond drops to zero abruptly and transfers the load to its nearest neighbor bonds.An explicit velocity-verlet algorithm was applied in solving the equation of motion of material points in the discrete PDLM model.e paper is organized as follows: we firstly reviewed the bond-based PD laminate theory in Section 2; then, the PD lattice model for 2D composite lamina was proposed in Section 3; the corresponding crack propagation experiment was introduced in Section 4; furthermore, the numerical validation scheme was implemented in a Fortran code, and benchmark problems including one elastic verification problem and two fracture problems of composite lamina were used to verify the model in Section 5; in the end of this paper, some conclusions were drawn based on the comparison of PD results with those of experiments or analytics.

Bond-Based Peridynamic Laminate Formulation
To remove limitations of continuum mechanics in dealing with discontinuities, peridynamic theory (PD) uses displacement integrations instead of differentials in its motion of equations.It assumes its reference configuration in a 3D Euclidean space, and at time t, the motion of equations is expressed as where f is the pairwise force scalar that material point x ′ (neighborhood of x ) exerts on material point x and H x represents the neighborhood of material point x.
For each material point x in the reference configuration (named P), its position changes from X to x under certain loads, and its neighbor material point P ′ deforms from X ′ to x ′ , as shown in Figure 1. e relative position vector ξ � x ′ − x is called a "bond," and the relative displacement vector is written as η � u ′ − u, in which u � x − X and u ′ � x ′ − X ′ are displacement vectors of material points P and P ′ , respectively.
To express the nonlocal interaction of a material point, a non-negative parameter δ, called horizon, is introduced to determine the interaction range of its family nodes, ( It is assumed that in bond-based PD, pairwise force depends only on ξ and ξ + η: where ‖ • ‖ denotes the 2-norm of a vector.ere exist two constraints on PD formulation; the first is that it must satisfy Newton's third law, And the second constraint comes from the law of conservation of angular momentum, 2 Advances in Materials Science and Engineering ( For the material point x ′ within the neighborhood of material point x, the magnitude of the pairwise force is expressed as where μ(t, ξ) is a history-dependent Boolean function which is either 0 (bond breakage) or 1(bond working).c is the material micromodulus, which can be de ned as a constant or a variable value that depends on the bond length, and it can be obtained by strain energy equivalence; s is the bond stretch, a notion that is similar to strain It should be noted that the de nition of bond stretch does not depend on time t or other history-relevant variables like fatigue and viscoplastic, and the material is regarded as nonmemory material; Silling [5] de nes such a microelastic material, the relationship between pairwise force and pairwise potential is expressed as e PD laminate model [6] describes normal and shear deformation of composite laminate with four types of bond: in-plane normal and shear bonds, and interlayer normal and shear bonds.For 2D composite lamina, only in-plane bonds are needed, as depicted in Figure 2. Most of the published PD composite unidirectional lamina models are bond based; in these models, because sti ness in ber direction is quite larger than that in the transverse direction, material anisotropic is characterized by distinguishing the bond properties from ber direction and others.Bonds parallel to ber direction are called ber bonds, and bonds in all other directions are called matrix bonds.
e bond micromodulus can be written in a uniform way as follows: where Δ denotes the Dirac delta function, m and f denote the matrix and ber, separately, and ξ 2 represents the ber direction vector.e constitutive model of composite lamina is illustrated in Figure 2; ber and matrix bonds deform elastically until their stretches reach to a critical value, then they fail abruptly and do not bear load anymore; s is the short of the stretch, t and c are the short of tension and compression.

Formulation of Bond-Based Peridynamic
Lattice Model for Composite Lamina e lattice model was adopted into bond-based PD framework for composite lamina to gain a higher computation e ciency by prede ning neighborhood of material points, while maintaining the capability of the PD model in presenting nonlinear behavior and damage.Another advantage of lattice methods is that they are sensitive to speci c directions which make them suitable for composite modeling.
eoretical derivation of the peridynamic lattice model (PDLM) is given at rst, and its discretized motion of equations is solved by an explicit time integration algorithm; then, the model was tested by the elastic response of a square plate and nally applied to the fracture prediction of composite lamina under tension in Section 5.
Various lattice forms can be chosen [23] for the PDLM, for example, triangular, rectangular, and honeycomb; in this paper, triangular lattice accounting for nearest neighbors is adopted as depicted in Figure 3, and angular forces between bonds are left out (bonds are of spring form); R is the radius of a lattice cell, n i , i 1, 2, . . ., 6, is the bond direction unit vectors.
It is assumed [11] that strain energy in the longitudinal direction is stored in ber bonds only, while strain energy along transverse is stored in matrix bonds; this assumption is feasible in continuous PD method; however, it is not applicable for PDLM.Given a homogeneous deformation s, the equivalence of strain energy (density) of the ber bonds (for example, ber bond directions n 1 and n 4 ) and the strain energy of homogenized orthotropic material can be obtained.And, the same is for the transverse direction.e strain energy density of the PDLM in ber and matrix directions can be expressed as Advances in Materials Science and Engineering where w is the micropotentials of each bond.To get an accuracy evaluation of bonds in PDLM, similar to the procedure by Askari et al. [15], a spherical function was applied to distinguish the matrix bond direction in a uniform way; the micropotentials can be written as where, c f and c m are the micromodulus of ber and matrix bond and D is a homogeneous factor dependent of ber bond direction vector (angle θ), which is chosen to be a Legendre function D(θ) ∞ n 0 A 0 n P 0 n (cos θ); the second and the fourth orders are selected, and e orthotropic property of composite is represented by rotating the lattice direction, which is a common method in lattice-based methods [24], as shown in Figure 4; ber bond direction is coincident with the local coordinate.
rough some mathematical manipulations, and assuming Poisson's ratio to be 1/3 (this limitation can be overcome by considering the volume potential energy of the cell, and will be discussed in another paper), micromodulus of a bond can be obtained, and detailed derivation process can be seen in, Fiber and matrix bond breakage is determined by a critical value s c ; s c is obtained by equating energy that breaks all the bonds across the fracture surface to the energy release rate separating halves of the body.Following the procedure of Silling and Askari [7], where G 1 ICr and G 2 ICr are energy release rates of fracture modes I and II .e criterion is expressed as Local damage of total, ber, and matrix can be expressed as e lattice method itself is an excellent numerical method for dynamic problems; it prede nes positions for integration (Gaussian) points of di erent subdomains, and then the integral form motion of equations of PDLM can be converted into a nite sum: where N is the number of cells(A lattice) in the body and N L ≤ 6 is the number of points in a lattice.An explicit leapfrog integration method, with secondorder accuracy, is applied in PDLM numerical simulations: where v i and f i /ρ are the velocity and acceleration vectors of material point i, Δt is the timestep, and v m+(1/2) i stands for the velocity at time t m + (1/2).

Experimental Research on Crack Propagation of Composite Lamina
In order to verify the validity and accuracy of the PDLM, tension tests of unidirectional (UD) composite plates with defects under tension were carried out.Failure modes of composite unidirectional plates with di erent defect (precrack) forms and di erent ber directions are obtained, which will provide a comparison and reference for the subsequent PDLM simulation.In this section, we rstly introduced the test equipment and samples and then introduced test process in detail and evaluated the test results qualitatively.

Test Equipment and Samples.
e test equipment used in the test is Zwick/Z100 electronic universal testing machine, as shown in Figure 5. e tension test follows the standard ASTM D5766.
e composite sample used in the experiment is a little di erent from that in the standard ASTM D5766.e total length of the sample is l 120 mm, e ective working length is w 70 mm, width is b 40 mm, and thickness is t 0.6 mm.Two forms of defects (precracks), i.e., hole and line cracks are considered.e schematic diagram and pictures of two samples are shown in Figure 6, Tension displacement loadings are applied at the clamped boundaries.According to the form of defect and the direction of tension of the unidirectional plate, the samples can be categorized into two types: open hole (A) and line crack (B).eir full names contain ber orientation α, and plate with line crack also contains a line angle β; the de nition of these two angles are shown in Figure 7, for example, B-45-30 indicates a plate with a hole, its ber direction is 45 °, and the line crack angle is 30 °.  e composite material is Torayca T700-12K (UD-T700-100gsm-38%), and speci c material parameters are shown in Table 1.Samples used in the experiment are listed in Table 2.  Advances in Materials Science and Engineering the loading rate of A0 and B0 specimens is set to 1 mm/min, and that of other types of specimens is 0.1 mm/min, which is uniformly loaded until the specimens are completely broken.

Test Process and Result Analysis.
Experimental results of three A-type specimens are given in Figure 8.At rst, for A-0 specimen, it can observed that the damage mechanism is ber breakage and matrix crack, and some irregular ber crack can be seen at the loading ends (red lines indicate the crack path); for A-45 specimen, the existence of precrack makes the sti ness of specimens change slightly during loading, and the ultimate failure is instantaneous brittle fracture.In this case, the strength of the sample is mainly determined by the interfacial properties between the bers and the matrix.
e failure mode is mainly the tensile cracking of the matrix, and generally there is no ber fracture, as shown in Figure 8(b).e nal fracture surfaces of the specimens are very at, and the cracks begin at the defect position of the specimens and propagate strictly along the direction of the bers.e rst place where matrix cracks occur is about 2 mm away from the top and bottom of the hole to the axis of the hole center, they are centrally symmetrical with respect to the hole center, and the cracks propagate in the opposite direction along the direction of the ber; for A-90, the brittle crack path is along the ber direction, and crack initiation position is about 1.1 mm from the vertical centerline of the hole, as shown in Figure 8(c).
e crack propagation results of B-type specimens with α 90 °are illustrated in Figure 9; the failure modes of all test specimens are matrix cracking.Cracks begin at the tip of the initial crack and propagate along the direction of the bers, leading to the ultimate failure and fracture.

Numerical Implementation and Validation
e PDLM can be still regarded to be a meshless particle method as bond-based PD [7], which discretizes the material into a series of arranged material points.ese material points contain all physical information of the material.Taking each material point as the object of study, the movement of material points at a certain time can be determined by calculating the interaction between material points, i.e., bonds.e motion state of all material points re ects the movement of the macromodel.In this section, we will show the code implementation process and then use this code to solve benchmark problems.

Code Implementation.
e code implementation process is illustrated in Figure 10, the di erence between bond-based PD and PDLM is that (1) the spatial discretion is xed before the time integration; (2) boundary and initial conditions are applied directly; and (3) when code enters into time integration loop, spatial integration does not need anymore; however, mass matrix should still be reassembled.After code implementation nishes, elastic response problems of composite lamina are applied to verify the accuracy of the model; then, two composite lamina plates with precracks under tension are solved, their results are compared to that of experiment, more attention is paid to the 45-degree crack propagation result, which shows the characteristics of lattice based models; at last, a convergence study of PDLM is conducted, and we pay main focus on the e ect of lattice size (radius) on the crack propagation results.

Elastic Response of Composite Unidirectional Reinforced
Lamina under Tension.After implementing the PDLM model in a Fortran code, the model was applied on the problem of a unidirectional reinforced lamina under tension.e fiber angle with respect to the horizontal coordinate is α, and three different α ′ s are considered, i.e., 0 ∘ , 45 ∘ , and 90 ∘ ; the discretized model by lattice at different angles can be seen in Figure 11.Lattice size is chosen to be R � 1 mm.Specimen geometry and boundary conditions setting is the same to that in Figure 7. Properties of the unidirectional reinforced lamina are shown in Table 1.e elastic response of the specimen without precracks is obtained to validate the accuracy of the PDLM simulation.
e horizontal displacement of the plate with different fiber angles is shown in Figure 12 (unit of legend bar is m).
To validate the accuracy of the model, we compare displacement results on the centerline nodes in the specimen with fiber angle 0 °with analytical results and present them in Figure 13.It is observed that displacement obtained by PDLM agrees well with analytical results.However, there is a little difference near the end of the plate due to the boundary enhancement deficiency of the meshfree method itself, and this deficiency can be alleviated by finer mesh.

Progressive Failure Process of Composite Unidirectional
Reinforced Lamina under Tension.After verification of the elastic response of undamaged plate, composite unidirectional reinforced lamina with two different precracks is considered, and attention will be paid on their fracture mode and progressive failure process.A constant velocity v � 2.0 × 10 −3 m/s was applied at the boundary nodes, as shown in Figure 6. e circle radius is r � 5 mm, and line crack length is 2a � 8 mm.As to the damage, according to several experimental results [25,26], transverse matrix cracking dominates the splitting crack under tension; thus only, matrix damage is considered.8

Advances in Materials Science and Engineering
Advances in Materials Science and Engineering en, the progressive failure process of A-0 specimen is given in Figure 14; matrix crack initiates at the top and bottom of precrack at about 20 μs and propagates towards ber direction, which corresponds to the ber splitting phenomenon as shown in Figure 8(a).As time increases, matrix crack nally propagates to the ends, the specimen is separated into three parts, which presents a ber/matrix separation failure pattern, and some microcracks were observed at the ends.e nal matrix damage, ber damage, and total damage are shown in Figure 15.
Comparison of PDLM results with experimental results is shown in Figure 16, and good agreement is observed.And, there is obvious ber damage in both results.Due to the randomness of experiment setup and specimen preparation, the fracture surface at ends is shown to be asymmetric, which should be symmetric analytically though.
Similarly, for the B-0-0 specimen, the progressive failure process predicted by PDLM is depicted in Figures 17 and 18. Crack also initiates at the end of the precrack and propagates in the ber direction.As time increases, ber damage emerges at the clamped ends.
Next, we investigate the damage pattern of A-45/90, B-45/90-90 specimens; the progressive failure process by PDLM is presented in Figure 19.Damage is dominated by Advances in Materials Science and Engineering          10 Advances in Materials Science and Engineering matrix damage, crack initiates at about 2 mm away from the vertical centerline of the hole in each A-type specimen and initiates at the precrack tip in each B-type specimen, and cracks nally propagate to the clamped ends, which show good agreement with experimental results shown in Figure 19.e agreement between PDLM and experimental results shows that PDLM can model complex damage forms, while pure lattice methods may fail to predict [19].

Progressive Process of Composite Unidirectional Reinforced Lamina under Mixed-Mode Loading.
In this part, we consider the failure pattern of unidirectional reinforced lamina under mixed-mode loading via change of line crack angle. is experiment was done on unidirectional reinforced lamina with ber direction α 90 °, and all other conditions are the same with above model.PD crack propagation results are compared with those of experiments as depicted in Figure 20.
Because ber strength is stronger than that of matrix, matrix damage dominates the total damage.As depicted in Figure 20, crack initiates at the ends of the precrack and propagates along the ber direction to boundaries, nally splitting the specimen into two parts, which coincides with experimental results very well.

Conclusions
In this paper, we proposed a hybrid peridynamics and lattice method (PDLM) formulation for simulating fracture in a UD composite lamina.Analytical derivations of material properties in PDLM for ber and matrix bond are obtained by a homogeneous assumption.Bond micromodulus is obtained by the equivalence of strain energy density to that in classical laminated theory, and critical stretch is related to modes I and II of fracture energies.Elastic response of composite lamina with ber directions of 0 ∘ , 45 ∘ , and 90 ∘ under tension is captured by PDLM.e results agree well with those of experiment.Numerical analysis is also performed for lamina with circle and line precracks, and crack propagation results show a good agreement with experimental observations.Progressive     Advances in Materials Science and Engineering failure process of lamina under mix-mode loading is also obtained by changing directions of the line crack.e PDLM model introduced in this paper focuses on 2D UD composite lamina; however, it can be easily extended to 3D composite laminates by adopting 3D lattice in the future to study fracture pattern of various ber orientations.Furthermore, the present PDLM model is still limited to quasi-isotropic lay-ups with 3 xed ber angles, and modeling composite with arbitrary ber direction will be developed in the future.

Appendix
In order to present the anisotropy of composites and take into account the distribution of bers and matrix, it is assumed that the micromodulus in the direction of bers has one value, while the other micromodulus from the direction of bers has another value, as shown in Figure 1.When the bers are in the positive direction along the axis, the bond's micromodulus can be written as where c f and c m are the ber and matrix bond micromoduli, respectively.As shown in Figure 21, composite lamina lies in the coordinate plane, and coordinate x denotes the ber direction.To characterize the continuous changes of ber bonds, a spherical function is introduced:   (A.8) According to the theory of composite mechanics, the stiffness coefficient of a laminate plate in any direction is as follows: It can be seen that the maximum number of trigonometric functions is four times.In order to establish the relationship with classical layer theory and better reflect the relationship between the micromodulus function and the angle, the expansion of c f (ϕ) is taken to the fourth order as follows: After simplification, we obtain the following form: (A.12) In order to obtain the relationship between micromodulus function and engineering constants such as elastic modulus, energy equivalence is applied.In the classical theory of elasticity, the stress-strain relationship of orthotropic material under plane stress can be written as Advances in Materials Science and Engineering where Q ij is the stiffness matrix coefficient and e strain energy density of an element under strain ε i is where, i, j � 1, 2, 6 and ε 6 � c 12 .In order to derive the relationship between the bond's modulus function and the stiffness matrix coefficient, four equations are constructed by assuming four different strain states and then solved simultaneously.e strain state and the corresponding bond elongation are shown in Table 3.
In order to obtain the bond stretch, the strain vector [ε 1 , ε 2 , c 12 ] is transformed into a new coordinate system: s � ε 1 cos 2 φ + ε 2 sin 2 φ + c 12 cos φ sin φ. (A.16) Under the condition of plane stress, the strain energy density of a material point is ξt dξ dφ.(A.17) In PD theory, the strain energy densities under four strain states are as follows: n � 1, With these coefficients obtained, consider the relation between elastic potential and bond length as follows [27]:

Figure 6 :Figure 7 :
Figure 6: Pictures and schematic of two plates with defects.(a) Plate with a hole.(b) Plate with a line crack.

Figure 10 :
Figure 10: Flow chart of the PDLM algorithm.

Table 1 :
Material parameters of composite lamina.

Table 2 :
Samples used in the experiment.

Table 3 :
Four strain states and the corresponding stretches.� ξ, other components are 0 s � ξ cos 2 φ 2 ε 2 � ξ, other components are 0 s � ξ sin 2 φ 3 c 12 � ξ, other components are 0 s � ξ cos φ sin φ 4 ε 1 � ε 2 � ξ, other components are 0 s � ξ 14 Advances in Materials Science and Engineering where E tot is the total energy density in a lattice cell and ξ i , i � 1, 2, ..., 6, is the length of the ith bond; then combining equations (A.18), (A.23), and (A.25), and assuming ϕ � 0 ∘ and 60 ∘ , we obtain the expression of micromodulus for a lamina: It can be seen that, when E 11 � E 22 , i.e., the material is isotropic, we can get c f � c m .Equation (A.26) degenerates to the original PD model for isotropic materials.