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.
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 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.
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 time t, the motion of equations is expressed aswhere is the pairwise force scalar that material point (neighborhood of ) exerts on material point and represents the neighborhood of material point .
For each material point in the reference configuration (named P), its position changes from to under certain loads, and its neighbor material point deforms from to , as shown in Figure 1. The relative position vector is called a “bond,” and the relative displacement vector is written as , in which and are displacement vectors of material points P and , 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. There 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,
For the material point within the neighborhood of material point , the magnitude of the pairwise force is expressed aswhere is a history-dependent Boolean function which is either 0 (bond breakage) or 1(bond working). 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; is the bond stretch, a notion that is similar to strain
It should be noted that the definition 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] defines such a microelastic material, the relationship between pairwise force and pairwise potential is expressed as
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:where denotes the Dirac delta function, and denote the matrix and fiber, separately, and 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; is the short of the stretch, and are the short of tension and compression.
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 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, , is the bond direction unit vectors.
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 , the equivalence of strain energy (density) of the fiber bonds (for example, fiber bond directions and ) 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 aswhere 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 aswhere, and are the micromodulus of fiber and matrix bond and is a homogeneous factor dependent of fiber bond direction vector (angle ), which is chosen to be a Legendre function ; the second and the fourth orders are selected, and . 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.
Through some mathematical manipulations, and assuming Poisson’s ratio to be (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 ; 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 and are energy release rates of fracture modes I and II . The criterion is expressed as
Local damage of total, fiber, and matrix can be expressed as
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:where is the number of cells(A lattice) in the body and is the number of points in a lattice.
An explicit leapfrog integration method, with second-order accuracy, is applied in PDLM numerical simulations:where and are the velocity and acceleration vectors of material point , is the timestep, and stands for the velocity at time .
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 Figure 5. The tension test follows the standard ASTM D5766.
The composite sample used in the experiment is a little different from that in the standard ASTM D5766. The total length of the sample is , effective working length is , width is , and thickness is . 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°.
(a)
(b)
The composite material is Torayca T700-12K (UD-T700-100gsm-38%), and specific material parameters are shown in Table 1. Samples used in the experiment are listed in Table 2.
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 Figure 8. 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).
(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.
(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 Figure 10, 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.
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 are considered, i.e., ; the discretized model by lattice at different angles can be seen in Figure 11. Lattice size is chosen to be . 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.
(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 Figure 12 (unit of legend bar is m).
(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 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.
(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 velocity was applied at the boundary nodes, as shown in Figure 6. The circle radius is , and line crack length is . 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 Figure 14; matrix crack initiates at the top and bottom of precrack at about 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.
(a)
(b)
(a)
(b)
(c)
Comparison of PDLM results with experimental results is shown in Figure 16, 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.
(a)
(b)
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 fiber direction. As time increases, fiber damage emerges at the clamped ends.
(a)
(b)
(c)
(d)
(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 Figure 19. 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].
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 , 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 fiber 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 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 of , , and 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.
Appendix
In order to present the anisotropy of composites and take into account the distribution of fibers and matrix, it is assumed that the micromodulus in the direction of fibers has one value, while the other micromodulus from the direction of fibers has another value, as shown in Figure 1. When the fibers are in the positive direction along the axis, the bond’s micromodulus can be written aswhere and are the fiber and matrix bond micromoduli, respectively. As shown in Figure 21, composite lamina lies in the coordinate plane, and coordinate x denotes the fiber direction.To characterize the continuous changes of fiber bonds, a spherical function is introduced:where is a jointed Legendre function,
(a)
(b)
Two-dimensional composite unidirectional plates can be regarded as transversely isotropic materials, i.e., the micromodulus function is independent of . Thus, when , we have , and then c has the following form:
At the same time, is symmetric along the plane; thus,
Using the properties of the jointed Legendre function
Obviously, when , . Then, can be expanded aswhere
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 is taken to the fourth order as follows:
Then,
After simplification, we obtain the following form:
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 aswhere is the stiffness matrix coefficient and
The strain energy density of an element under strain iswhere, and . 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. The strain state and the corresponding bond elongation are shown in Table 3.
In order to obtain the bond stretch, the strain vector is transformed into a new coordinate system:
Under the condition of plane stress, the strain energy density of a material point is
In PD theory, the strain energy densities under four strain states are as follows:
n = 1,
n = 2,
n = 3,
n = 4,
In conjunction with the strain energy density calculated by classical elasticity theory, the coefficients can be expressed as
Finally, the micromodulus is expressed aswhere
With these coefficients obtained, consider the relation between elastic potential and bond length as follows [27]:where is the total energy density in a lattice cell and , is the length of the ith bond; then combining equations (A.18), (A.23), and (A.25), and assuming and , we obtain the expression of micromodulus for a lamina:
It can be seen that, when , i.e., the material is isotropic, we can get . Equation (A.26) degenerates to the original PD model for isotropic materials.
Data Availability
The numerical and experimental data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.
Acknowledgments
This research was supported by the Natural Science Foundation of China (Grant no. 11772110). And, the original Fortran codes come from Madenci and Oterkus (2014). These supports are gratefully acknowledged.