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

Authors: Jushang Guo; Weicheng Gao; Zhenyu Liu; Xiongwu Yang; Fengshou Li
Journal: Advances in Materials Science and Engineering
 (2019)

---

## Abstract

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. The peridynamic lattice model is found to be robust and successful in modeling dynamic brittle fracture of 2D composite lamina and can be extended to composite laminates by applying 3D lattice structure.


---

## Body

## 1. 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. Then, computation savings can be obtained by predefining neighbor of particles in a fixed topology before entering time integration loop. The 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–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. The 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.The paper is organized as follows: we firstly reviewed the bond-based PD laminate theory in Section2; 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.


## 2. 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 timet, the motion of equations is expressed as(1)ρu¨=∫Hxfη,ξdVx′+bx,t,where f is the pairwise force scalar that material point x′ (neighborhood of x) exerts on material point x and Hx represents the neighborhood of material point x.For each material pointx 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. The 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.Figure 1
Schematic of motion of peridynamic material points.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,(2)fη,ξ=0,ifξ>δ.It is assumed that in bond-based PD, pairwise force depends only onξ and ξ+η:(3)fη,ξ=fη,ξη+ξη+ξ,where ⋅ denotes the 2-norm of a vector. There exist two constraints on PD formulation; the first is that it must satisfy Newton’s third law,(4)fη,ξ=−f−η,−ξ.And the second constraint comes from the law of conservation of angular momentum,(5)η+ξ×fη,ξ=0,∀η,ξ.For the material pointx′ within the neighborhood of material point x, the magnitude of the pairwise force is expressed as(6)f=μt,ξcsη,ξ,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 defined 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(7)sη,ξ=η+ξ−ξξ.It should be noted that the definition of bond stretch does not depend on timet or other history-relevant variables like fatigue and viscoplastic, and the material is regarded as nonmemory material; Silling [5] defines such a microelastic material, the relationship between pairwise force and pairwise potential is expressed as(8)fη,ξ=∂wη,ξ∂η.The 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 stiffness in fiber direction is quite larger than that in the transverse direction, material anisotropic is characterized by distinguishing the bond properties from fiber direction and others. Bonds parallel to fiber direction are called fiber bonds, and bonds in all other directions are called matrix bonds. The bond micromodulus can be written in a uniform way as follows:(9)cξ=cmξ+cfξΔξ2,where Δ denotes the Dirac delta function, m and f denote the matrix and fiber, separately, and ξ2 represents the fiber direction vector. The constitutive model of composite lamina is illustrated in Figure 2; fiber 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.Figure 2
Constitutive model of composite lamina.


## 3. Formulation of Bond-Based Peridynamic Lattice Model for Composite Lamina

The lattice model was adopted into bond-based PD framework for composite lamina to gain a higher computation efficiency by predefining 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 specific directions which make them suitable for composite modeling. Theoretical derivation of the peridynamic lattice model (PDLM) is given at first, 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 finally applied to the fracture prediction of composite lamina under tension in Section5.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, ni,i=1,2,…,6, is the bond direction unit vectors.Figure 3
Triangular lattice accounting for nearest neighbors. (a) Cell representation. (b) Normal vector representation.It is assumed [11] that strain energy in the longitudinal direction is stored in fiber 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 fiber bonds (for example, fiber bond directions n1 and n4) and the strain energy of homogenized orthotropic material can be obtained. And, the same is for the transverse direction. The strain energy density of the PDLM in fiber and matrix directions can be expressed as(10)W1=Wf=12∫Hxwfn1+wfn4dAξ=12∫−llwfdAξ,W2=Wm=12∫HxwmnidAξ,i=2,3,5,6,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(11)wf=12cfs2ξ,wm=12Dθcms2ξ,where, cf and cm are the micromodulus of fiber and matrix bond and D is a homogeneous factor dependent of fiber bond direction vector (angle θ), which is chosen to be a Legendre function Dθ=∑n=0∞An0Pn0cosθ; the second and the fourth orders are selected, P20cosθ=1/23cos2θ−1 and P40cosθ=1/835cos4θ−30cos2θ+3. The 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; fiber bond direction is coincident with the local coordinate.Figure 4
Fiber direction represented by lattice rotation.Through some mathematical manipulations, and assuming Poisson’s ratio to be1/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,(12)cf=14.34E11−8.34E22πtR3,cm=6.09E11−0.09E22πtR3.Fiber and matrix bond breakage is determined by a critical valuesc; sc 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],(13)scf=4πGICr19E1R,scm=4πGICr29E2R,where GICr1 and GICr2 are energy release rates of fracture modes I and II . The criterion is expressed as(14)μf=0,sf≥scf,1,others,μm=0,sm≥scm,1,others.Local damage of total, fiber, and matrix can be expressed as(15)ϕtot=1−∫HxφdA∫Hx1⋅dA,ϕf=1−∫Hx∩HfφfdA∫Hx∩Hf1⋅dA,ϕm=1−∫Hx∩HmφmdA∫Hx∩Hm1⋅dA.The lattice method itself is an excellent numerical method for dynamic problems; it predefines positions for integration (Gaussian) points of different subdomains, and then the integral form motion of equations of PDLM can be converted into a finite sum:(16)ρ0y¨b=∑i=1N∑j=1NLfη,ξ,t+bb,where N is the number of cells(A lattice) in the body and NL≤6 is the number of points in a lattice.An explicit leapfrog integration method, with second-order accuracy, is applied in PDLM numerical simulations:(17)vim+1/2=vim+Δt2ρifim,uim+1=uim+vim+1/2Δt,vim+1=vim+1/2+Δt2ρifim+1,where vi and fi/ρ are the velocity and acceleration vectors of material point i, Δt is the timestep, and vim+1/2 stands for the velocity at time t=m+1/2.


## 4. 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 different defect (precrack) forms and different fiber directions are obtained, which will provide a comparison and reference for the subsequent PDLM simulation. In this section, we firstly introduced the test equipment and samples and then introduced test process in detail and evaluated the test results qualitatively.


### 4.1. Test Equipment and Samples

The test equipment used in the test is Zwick/Z100 electronic universal testing machine, as shown in Figure5. The tension test follows the standard ASTM D5766.Figure 5
The Zwick/Z100 test machine.The composite sample used in the experiment is a little different from that in the standard ASTM D5766. The total length of the sample isl=120mm, effective working length is w=70mm, width is b=40mm, and thickness is t=0.6mm. Two forms of defects (precracks), i.e., hole and line cracks are considered. The 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). Their full names contain fiber orientation α, and plate with line crack also contains a line angle β; the definition of these two angles are shown in Figure 7, for example, B-45-30 indicates a plate with a hole, its fiber direction is 45°, and the line crack angle is 30°.Figure 6
Pictures and schematic of two plates with defects. (a) Plate with a hole. (b) Plate with a line crack.
(a)
(b)Figure 7
Definition of the fiber direction angle and the line crack angle.The composite material is Torayca T700-12K (UD-T700-100gsm-38%), and specific material parameters are shown in Table1. Samples used in the experiment are listed in Table 2.Table 1
Material parameters of composite lamina.
Young’s modulus (MPa)
Poisson’s ratio
Shear modulus (MPa)
E
1
E
2
E
3
ν
12
ν
13
ν
23
G
12
G
13
G
23
105000
8400
8400
0.32
0.3
0.3
4000
3500
3500Table 2
Samples used in the experiment.
A
A-0
A-45
A-90
B
B-0-0/30/45/60
B-45-0/30/45/60
B-00-0/30/45/60


### 4.2. Test Process and Result Analysis

All specimens are loaded by displacement loadings. Because the unidirectional strength of fibers is much greater than the vertical strength of fibers, 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.Experimental results of three A-type specimens are given in Figure8. At first, for A-0 specimen, it can observed that the damage mechanism is fiber breakage and matrix crack, and some irregular fiber crack can be seen at the loading ends (red lines indicate the crack path); for A-45 specimen, the existence of precrack makes the stiffness 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 fibers and the matrix. The failure mode is mainly the tensile cracking of the matrix, and generally there is no fiber fracture, as shown in Figure 8(b). The final fracture surfaces of the specimens are very flat, and the cracks begin at the defect position of the specimens and propagate strictly along the direction of the fibers. The first 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 fiber; for A-90, the brittle crack path is along the fiber direction, and crack initiation position is about 1.1 mm from the vertical centerline of the hole, as shown in Figure 8(c).Figure 8
Final failure morphology of three A-type composite specimens. (a)α = 0°. (b) α = 45°. (c) α = 90°.
(a)
(b)
(c)The 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 fibers, leading to the ultimate failure and fracture.Figure 9
Final failure morphology of four B-type composite specimens. (a)β=0°. (b) β=30°. (c) β=45°. (d) β=60°.
(a)
(b)
(c)
(d)


## 4.1. Test Equipment and Samples

The test equipment used in the test is Zwick/Z100 electronic universal testing machine, as shown in Figure5. The tension test follows the standard ASTM D5766.Figure 5
The Zwick/Z100 test machine.The composite sample used in the experiment is a little different from that in the standard ASTM D5766. The total length of the sample isl=120mm, effective working length is w=70mm, width is b=40mm, and thickness is t=0.6mm. Two forms of defects (precracks), i.e., hole and line cracks are considered. The 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). Their full names contain fiber orientation α, and plate with line crack also contains a line angle β; the definition of these two angles are shown in Figure 7, for example, B-45-30 indicates a plate with a hole, its fiber direction is 45°, and the line crack angle is 30°.Figure 6
Pictures and schematic of two plates with defects. (a) Plate with a hole. (b) Plate with a line crack.
(a)
(b)Figure 7
Definition of the fiber direction angle and the line crack angle.The composite material is Torayca T700-12K (UD-T700-100gsm-38%), and specific material parameters are shown in Table1. Samples used in the experiment are listed in Table 2.Table 1
Material parameters of composite lamina.
Young’s modulus (MPa)
Poisson’s ratio
Shear modulus (MPa)
E
1
E
2
E
3
ν
12
ν
13
ν
23
G
12
G
13
G
23
105000
8400
8400
0.32
0.3
0.3
4000
3500
3500Table 2
Samples used in the experiment.
A
A-0
A-45
A-90
B
B-0-0/30/45/60
B-45-0/30/45/60
B-00-0/30/45/60


## 4.2. Test Process and Result Analysis

All specimens are loaded by displacement loadings. Because the unidirectional strength of fibers is much greater than the vertical strength of fibers, 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.Experimental results of three A-type specimens are given in Figure8. At first, for A-0 specimen, it can observed that the damage mechanism is fiber breakage and matrix crack, and some irregular fiber crack can be seen at the loading ends (red lines indicate the crack path); for A-45 specimen, the existence of precrack makes the stiffness 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 fibers and the matrix. The failure mode is mainly the tensile cracking of the matrix, and generally there is no fiber fracture, as shown in Figure 8(b). The final fracture surfaces of the specimens are very flat, and the cracks begin at the defect position of the specimens and propagate strictly along the direction of the fibers. The first 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 fiber; for A-90, the brittle crack path is along the fiber direction, and crack initiation position is about 1.1 mm from the vertical centerline of the hole, as shown in Figure 8(c).Figure 8
Final failure morphology of three A-type composite specimens. (a)α = 0°. (b) α = 45°. (c) α = 90°.
(a)
(b)
(c)The 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 fibers, leading to the ultimate failure and fracture.Figure 9
Final failure morphology of four B-type composite specimens. (a)β=0°. (b) β=30°. (c) β=45°. (d) β=60°.
(a)
(b)
(c)
(d)


## 5. Numerical Implementation and Validation

The 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. These 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. The motion state of all material points reflects the movement of the macromodel. In this section, we will show the code implementation process and then use this code to solve benchmark problems.


### 5.1. Code Implementation

The code implementation process is illustrated in Figure10, the difference between bond-based PD and PDLM is that (1) the spatial discretion is fixed 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 finishes, 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 effect of lattice size (radius) on the crack propagation results.Figure 10
Flow chart of the PDLM algorithm.


### 5.2. 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. The fiber angle with respect to the horizontal coordinate isα, and three different α′s are considered, i.e., 0∘,45∘, and90∘; the discretized model by lattice at different angles can be seen in Figure 11. Lattice size is chosen to be R=1mm. 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.Figure 11
Lattice mesh with different fiber angles. (a)α=0∘. (b) α=45∘. (c) α=90∘.
(a)
(b)
(c)The elastic response of the specimen without precracks is obtained to validate the accuracy of the PDLM simulation. The horizontal displacement of the plate with different fiber angles is shown in Figure12 (unit of legend bar is m).Figure 12
Horizontal displacement of UD reinforced lamina with different fiber directions. (a)α=0∘. (b) α=45∘. (c) α=90∘.
(a)
(b)
(c)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 Figure13. 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.Figure 13
Comparison of displacement between PDLM and analytical results. (a) Horizontal displacement. (b) Vertical displacement.
(a)
(b)


### 5.3. 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 velocityv=2.0×10−3m/s was applied at the boundary nodes, as shown in Figure 6. The circle radius is r=5mm, and line crack length is 2a=8mm. 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.Then, the progressive failure process of A-0 specimen is given in Figure14; matrix crack initiates at the top and bottom of precrack at about 20μs and propagates towards fiber direction, which corresponds to the fiber splitting phenomenon as shown in Figure 8(a). As time increases, matrix crack finally propagates to the ends, the specimen is separated into three parts, which presents a fiber/matrix separation failure pattern, and some microcracks were observed at the ends. The final matrix damage, fiber damage, and total damage are shown in Figure 15.Figure 14
Progressive failure process of A-0 specimen with a hole under tension: (a) displacement (m) and (b) damage (0-1).
(a)
(b)Figure 15
Final failure morphology of A-0 specimen with a hole under tension (damage (0-1)). (a) Matrix damage. (b) Fiber damage. (c) Total damage.
(a)
(b)
(c)Comparison of PDLM results with experimental results is shown in Figure16, and good agreement is observed. And, there is obvious fiber 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.Figure 16
Failure morphology of UD reinforced lamina with circle crack. (a) PD result. (b) Experiment.
(a)
(b)Similarly, for the B-0-0 specimen, the progressive failure process predicted by PDLM is depicted in Figures17 and 18. Crack also initiates at the end of the precrack and propagates in the fiber direction. As time increases, fiber damage emerges at the clamped ends.Figure 17
Progressive failure process of B-0-0 specimen (damage (0-1)). (a)t=30μs. (b) t=50μs. (c) t=70μs. (d) t=100μs.
(a)
(b)
(c)
(d)Figure 18
Failure morphology of B-0-0 specimen (damage (0-1)). (a) Matrix damage. (b) Fiber damage. (c) Total damage.
(a)
(b)
(c)Next, we investigate the damage pattern of A-45/90, B-45/90-90 specimens; the progressive failure process by PDLM is presented in Figure19. Damage is dominated by 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 finally propagate to the clamped ends, which show good agreement with experimental results shown in Figure 19. The agreement between PDLM and experimental results shows that PDLM can model complex damage forms, while pure lattice methods may fail to predict [19].Figure 19
Progressive damage of lamina with fiber directionα = 45° and α = 90° (damage (0-1)).


### 5.4. 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. This experiment was done on unidirectional reinforced lamina with fiber 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.Figure 20
Crack propagation of UD reinforced lamina under mixed-mode loading (damage (0-1)).Because fiber strength is stronger than that of matrix, matrix damage dominates the total damage. As depicted in Figure20, crack initiates at the ends of the precrack and propagates along the fiber direction to boundaries, finally splitting the specimen into two parts, which coincides with experimental results very well.


## 5.1. Code Implementation

The code implementation process is illustrated in Figure10, the difference between bond-based PD and PDLM is that (1) the spatial discretion is fixed 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 finishes, 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 effect of lattice size (radius) on the crack propagation results.Figure 10
Flow chart of the PDLM algorithm.


## 5.2. 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. The fiber angle with respect to the horizontal coordinate isα, and three different α′s are considered, i.e., 0∘,45∘, and90∘; the discretized model by lattice at different angles can be seen in Figure 11. Lattice size is chosen to be R=1mm. 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.Figure 11
Lattice mesh with different fiber angles. (a)α=0∘. (b) α=45∘. (c) α=90∘.
(a)
(b)
(c)The elastic response of the specimen without precracks is obtained to validate the accuracy of the PDLM simulation. The horizontal displacement of the plate with different fiber angles is shown in Figure12 (unit of legend bar is m).Figure 12
Horizontal displacement of UD reinforced lamina with different fiber directions. (a)α=0∘. (b) α=45∘. (c) α=90∘.
(a)
(b)
(c)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 Figure13. 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.Figure 13
Comparison of displacement between PDLM and analytical results. (a) Horizontal displacement. (b) Vertical displacement.
(a)
(b)


## 5.3. 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 velocityv=2.0×10−3m/s was applied at the boundary nodes, as shown in Figure 6. The circle radius is r=5mm, and line crack length is 2a=8mm. 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.Then, the progressive failure process of A-0 specimen is given in Figure14; matrix crack initiates at the top and bottom of precrack at about 20μs and propagates towards fiber direction, which corresponds to the fiber splitting phenomenon as shown in Figure 8(a). As time increases, matrix crack finally propagates to the ends, the specimen is separated into three parts, which presents a fiber/matrix separation failure pattern, and some microcracks were observed at the ends. The final matrix damage, fiber damage, and total damage are shown in Figure 15.Figure 14
Progressive failure process of A-0 specimen with a hole under tension: (a) displacement (m) and (b) damage (0-1).
(a)
(b)Figure 15
Final failure morphology of A-0 specimen with a hole under tension (damage (0-1)). (a) Matrix damage. (b) Fiber damage. (c) Total damage.
(a)
(b)
(c)Comparison of PDLM results with experimental results is shown in Figure16, and good agreement is observed. And, there is obvious fiber 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.Figure 16
Failure morphology of UD reinforced lamina with circle crack. (a) PD result. (b) Experiment.
(a)
(b)Similarly, for the B-0-0 specimen, the progressive failure process predicted by PDLM is depicted in Figures17 and 18. Crack also initiates at the end of the precrack and propagates in the fiber direction. As time increases, fiber damage emerges at the clamped ends.Figure 17
Progressive failure process of B-0-0 specimen (damage (0-1)). (a)t=30μs. (b) t=50μs. (c) t=70μs. (d) t=100μs.
(a)
(b)
(c)
(d)Figure 18
Failure morphology of B-0-0 specimen (damage (0-1)). (a) Matrix damage. (b) Fiber damage. (c) Total damage.
(a)
(b)
(c)Next, we investigate the damage pattern of A-45/90, B-45/90-90 specimens; the progressive failure process by PDLM is presented in Figure19. Damage is dominated by 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 finally propagate to the clamped ends, which show good agreement with experimental results shown in Figure 19. The agreement between PDLM and experimental results shows that PDLM can model complex damage forms, while pure lattice methods may fail to predict [19].Figure 19
Progressive damage of lamina with fiber directionα = 45° and α = 90° (damage (0-1)).


## 5.4. 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. This experiment was done on unidirectional reinforced lamina with fiber 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.Figure 20
Crack propagation of UD reinforced lamina under mixed-mode loading (damage (0-1)).Because fiber strength is stronger than that of matrix, matrix damage dominates the total damage. As depicted in Figure20, crack initiates at the ends of the precrack and propagates along the fiber direction to boundaries, finally splitting the specimen into two parts, which coincides with experimental results very well.


## 6. 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 fiber 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 fiber directions of0∘, 45∘, and 90∘ under tension is captured by PDLM. The 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 failure process of lamina under mix-mode loading is also obtained by changing directions of the line crack. The 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 fiber orientations. Furthermore, the present PDLM model is still limited to quasi-isotropic lay-ups with 3 fixed fiber angles, and modeling composite with arbitrary fiber direction will be developed in the future.



---

*Source: 3748795-2019-02-04.xml*