Progressive failure analysis of CFRP composite laminates under uniaxial tension using a discrete element method

Abstract

This study presents a 3D Discrete Element Method (DEM) model for the progressive failure analysis of Carbon Fiber Reinforced Polymer (CFRP) composite materials subjected to uniaxial tensile loading. Particles in the model are packed and bonded in regular patterns (hexagonal or square). The relationship between the bond stiffness and material properties is established based on the average strain energy method. The random distribution of bond strengths calibrated from experiments with a variation of 30% and 10% following a normal distribution law is assigned to the bonds in 0∘ and 90∘ plies to capture random cracks, respectively. Tsai-Hill failure criterion is utilized for the calibration of bond strength of $\pm 45°$ plies to predict their failures in composite laminates. Quantitative and qualitative analyses were conducted for predicting the damage initiation and propagation of the cross-ply and Quasi-Isotropic (QI) composite laminates under tensile loading, respectively. Two interface stiffnesses were utilized in the failure prediction of cross-ply composite laminates, and it was found that the numerical results with the interface stiffness calculated from fracture energy are in good, quantitative agreements with the experiments. All the four stages of the failure process of QI composite laminates are well captured by the 3D DEM model, including isolated cracks, inner delamination cracks, outer delamination cracks and final failure.

Introduction

Glass or carbon fiber reinforced polymer (GFRP or CFRP) composite has been widely used in aerospace, mechanical and civil engineering for several decades, mainly due to its high stiffness-weight ratio. Besides, with proper design and optimization of the layer-up of several individual laminae, specific desired modulus and strength in different directions of the laminates can be achieved. However, due to the complexity of the physical microstructure of fiber-reinforced composite materials, the onset of damage does not cause the instantaneous failure of the entire structure, and there exists a progressive process from the damage initiation to the final structural collapse. ¹ The local failure such as fiber breakage, matrix cracks, fiber-matrix debonding, and interfacial delamination, may cause some permanent loss of integrity within composite laminates, resulting in the loss of the stiffness and strength of the laminates. Therefore, understanding of the failure mechanism, prediction of the initiation and progression of the damage as well as the final maximum loads the structure can withstand is of significant importance for assessing the performance of the composite laminates and optimization of the reliable design. ²

Transverse cracking and delamination are typical damages before the final catastrophic failure when the laminate is subjected to tensile loading. Transverse cracking can be initiated from the matrix cracking and/or fiber/matrix interfacial debonding caused by the high tensile stress in 90∘ ply. Delamination is usually initiated by the shear stress concentrating between the two adjacent plies due to the different material properties of these plies, under tensile loadings or free-edge stresses. These damage processes have been studied by the researchers for many decades, theoretically, experimentally, and numerically.

Pagano et al. ³ developed an analytical model to predict the stress of the cross-ply laminates based on the assumption of steady-state cracking, subjected to the thermal and thermo-mechanical loading conditions. Meanwhile, they found that the steady-state cracking has occurred since the full-width cracks have developed. Gudmundson and Zang ⁴ proposed a new model to predict the thermoelastic properties of the composite laminates containing plies with transverse matrix cracks. However, this model only considered the average stress-strain relationship at given matrix crack densities, and it cannot deal with the local stress distribution at each ply. Additionally, this model did not consider the crack initiation and propagation, making the model less flexible for the different dynamic loading conditions. Kobayashi et al.^5,6 conducted an analysis on the prediction of the formation of transverse crack evolution with energy release rate and average stress criteria based on Gudmundson and Zang ⁴ in quasi-isotropic laminates. It was found that the thickness of 90∘ plies and the modulus of adjacent plies have a significant influence on transverse cracking behavior. Okabe et al. ⁷ formulated an analytical model as a function of crack density based on the continuum damage mechanics for predicting the stiffness reduction of composite laminates. It was proved that the model is capable of predicting the stiffness of laminates containing damage caused by transverse cracks (or surface crack) with only mechanical properties of lamina and the lay-up configuration, and calculating the local stress distribution within a ply including transverse cracks as a function of the crack density. Wang et al. ⁸ proposed a model that consists of an application of fracture mechanics and a rational representation of material flaw distribution, which is a new concept of effective flaws to replace the conventional constant ply strength criterion. Onodera et al. ⁹ proposed a refined stress field model based on Okabe et al. ⁷ by considering the residual thermal strain for the prediction of the transverse crack density evolution. Ogi ¹⁰ proposed a theoretical model for predicting the transverse crack propagation within the cross-ply laminates based on a statistical distribution of microsized flaws, with an assumption that the initiation of transverse cracks only occurs through the ply thickness.

Several numerical methods have been proposed in the last decade for predicting and analyzing cracking in composites such as the finite element method (FEM)^11–18 and Extended FEM (X-FEM).^19–22 Within the framework of FEM, the stress distribution that is calculated to meet the failure criteria of a composite laminate is obtained by assuming that the material is not damaged, and the stress redistribution is calculated by the strain softening induced by the local damage, which is not considered in most criteria. X-FEM is very effective when dealing with cracks or even multiple cracks in homogenous isotropic materials, for which the enrichment functions are available. However, for complex composite material systems which have heterogeneous characteristics, the X-FEM can only work for some special case such as a delamination crack at a symmetric plane. ²³ Besides, it is difficult for these methods to consider the interaction between different damage modes or the coupling among multiple cracks. A promising alternative numerical method called the phantom node method was proposed by Hansbo and Hansbo ²⁴ to tackle the arbitrary cracking problem. They found that the arbitrary discontinuity in an element can be dealt with by adding an extra element on the top of that element, so that each of them can account for the stiffness and force contribution from its bisected physical domain, in which two domains are connected by linear or nonlinear springs, ²⁴ or cohesive failure tractions. ²⁵ The main advantage of this method is that the standard finite element shape function can be utilized within the framework of the phantom node method, rendering it completely compatible with standard FEM programs. However, adding extra elements changes the numerical model dimension dynamically, thus requiring full access to the source code, which is not always possible for commercial software. ²⁶ In order to ease this limitation, Ling et al. ²⁷ proposed an element, called augmented FEM (A-FEM), with double nodes dealing with the arbitrary intra-element discontinuity. Afterward, the A-FEM was successfully applied to predict the damage initiation and progression and failure prediction in the composite materials with an implementation of a subroutine in the commercial FE code ABAQUS/Explicit. ²⁸

Apart from the FEM modeling, the discrete element method (DEM) has been introduced for modeling the damage propagation in the composite using a bond-based approach. Crack propagation in composite laminates and stress-strain curve were simulated under the transverse tensile loading by DEM in literature,^29–32 which suggests the capability of the DEM model for predicting the damage initiation, propagation until the final failure strength by providing more details such as crack density and stiffness reduction, etc. However, for a test specimen subjected to tensile loading, transverse cracks are usually initiated from the free edges of the specimen, then propagate almost instantaneously through the full width of the specimen in the case of monotonic loading. ³³ Thus, two dimensional DEM models are not capable of capturing the damage propagation in width direction, and three-dimensional models is more attractive when it comes to the damage propagation in any direction. Recently, a 3D DEM model was proposed by Maheo et al. ³⁴ to demonstrate its ability to capture the damage propagation in composite without any preset paths. This model was then improved by implementing a cohesive contact model for modeling the interfacial delamination of composite materials. ³⁵ However, it is still difficult to model the angle-ply damage propagation with their current 3D DEM model due to its geometry limitation.

In the authors’ previous work, DEM models have been developed to study the progressive failure in single fiber composite, ³⁶ FRP laminae ²⁹ as well as FRP cross-ply laminates. ³² In their latest research, ³⁷ the authors were able to validate a 3D DEM model by the classical lamination theory and FEM for its elasticity, showing its capability of modeling the lamina, cross-ply and angle-ply laminates. This study expands upon the previous work by extending the developed 3D DEM model to analyze the damage initiation and progression in cross-ply laminates and QI composite laminates with an open hole. Firstly, the 3D DEM model was built up by combining different particle arrangement techniques for separate angle plies. The strength of the parallel bonds in different plies was calibrated from experiments and the Tsai-Hill failure criterion. Then case studies were performed to demonstrate the progression of transverse cracking and delamination in cross-ply and QI CFRP composite laminates with comparisons to experimental measurements and observations.

3D DEM modeling methodology

As reported in our previous work, ³⁷ the extended 2D DEM model is capable of modeling the elasticity of the general anisotropic materials. The detailed theory of DEM, parallel bond model and particle packing strategies, used in the present 3D DEM models, can be found in literature.^37,38 In contrast to previous studies, the damage process of a composite laminate is taken into consideration by means of the breakage of bonds. For simplicity and convenience, both of the bonds within a lamina and at the interface between different layers are modeled by the linear parallel bonds. The constitutive behavior of the parallel bond at contact is illustrated in Figure 1. It can be seen in Figure 1 that the mechanical behavior of the parallel bond is linearly elastic up to its tensile and shear strengths (σ_max and τ_max). Beyond which, bond breakage occurs, resulting in a sharp decrease of the stiffness and, consequently, the force. The maximum normal stress and shear stress can be determined by ${\sigma }_n=\frac{\overline{F_n}}{\overline{A}}$ and $\tau =\frac{\overline{F_s}}{\overline{A}}$ ( $\overline{A}$ is the cross-section area of particles), respectively. It should be noted that from the perspective of numerical modeling, one single composite lamina has the same mechanical performance as multiple-layer composite laminae under tensile loading. Therefore, for the sake of computational cost, one single layer is utilized to calibrate its bond strength according to the experimental data of composite laminae under tensile loading, such as 0∘ and 90∘ ply lamina. In addition, two particles are bonded together with a parallel bond, and a disk is formed with its plane perpendicular to the bond length direction when the bond breaks. The red disk can be found in Figure 1, while in the two case studies, different colors are selected for better distinguishing the different types of cracks. Taking advantage of this characteristic, the disks can be collected to describe the real-time cracks within the composite during the damage propagation. Of note, more sophisticated non-linear constitutive models such as the cohesive zone model and multiple-bond model could be adopted for the bonds, but that is out of the scope of the current study.

Figure 1.

Constitutive behavior of the parallel bond.

In DEM models, particles can bond together at contacts and separate when the bond strength or energy is exceeded. Thus, it is capable of simulating the motion of individual particles and the behavior of bulk material assembled randomly or regularly with particles through bonds at contacts using specific constitutive laws. Different from the traditional solutions, which are based on the strain and stress concepts, contact properties are the predominant parameters in a DEM model, together with the geometry parameters of the particles. Subject to external loading, when the strength or the fracture energy of a bond between particles is exceeded, the bond breaks. Hence, damage modes and their interactions emanate during the process of debonding of particles. The way that DEM discretizes the material domain provides an alternative to tackle the discontinuity problems, such as crack tip singularities and crack formulation criteria in the damage/failure prediction of composite materials, thanks to its naturally discontinuous representation of composite materials using assembling particles in different strategies. ³²

Modeling 0∘-, $30°$ -, $60°$ - and 90∘-ply lamina

In the extended 2D DEM model, a lamina is modeled as a unidirectionally reinforced composite lamina by assuming a plane stress state with a six-spring hexagonal packing as a base unit, shown in Figure 2. Figure 2(a) and (b) represent the base unit for 0∘, 60∘ plies and 30∘, 90∘ plies, respectively. Noting that due to the limitation of the hexagonal geometry, the base unit of 90∘ and 30∘ plies can be obtained from a rotation of 90 degree of 0∘ and 60∘ base units to make sure there are bonds locating at the targeted angle plies. The constructed lamina from these base units can be found in Figure 3, in which the particles are hidden except at the boundaries, and red and blue bonds represent the fiber and matrix bonds, respectively. Here the parallel bonds aligned with the fiber direction are called fiber bonds and all others bonds matrix bonds. It is important to emphasize that the fiber bond elastic properties are not matched to the actual fiber properties but are determined by the average strain energy method based on the experimental elastic properties of the composite lamina along the fiber direction. The same is true for the matrix bonds. The strengths of these bonds are calibrated based on the experimental tensile strengths of 0∘ and 90∘ ply laminae.

Figure 2.

Basic unit of the hexagonal packing for the composite lamina. (For interpretation the diameter of particles is reduced to half of their original ones.).

Figure 3.

Extended 2D hexagonal DEM model of composite lamina. (Particles are hidden except the ones of left and right boundaries with different colours; bonds with different colors represent different stiffnesses.). ³⁷ (a) 0 ply, (b) 30 ply, (c) 60 ply and (d) 90 ply.

Each contact between the two particles is modeled with a linear parallel bond, represented by a set of elastic springs with constant normal stiffness k_n and shear stiffness k_s, respectively. For the sake of better visualization of the relation between particles and bonds, the diameter of particles is reduced to half of their original ones. The diameter of particles can be determined by dividing the thickness of laminates by the number of layers, while the numbers of particles in length and width directions can be calculated by dividing the length and width by the particle diameter, respectively. The relationship between the spring constants and elastic constants is established according to the average strain energy method. ³⁹ Therefore, these spring constants can be related to the macro-scale elastic properties of the material, such as Young’s Modulus and Poisson’s ratio. In this approach, a general formula, which is based on the unit of Figure 2(a) for 0∘ and 60∘ plies, was adopted from the Liu and Liu ³⁹ with average strain energy method for anisotropic materials, as follows:

$$\{\begin{array}{l}{k}_{n1}=\frac{\sqrt{3}}{6}(3C_{11}+2C_{12}-C_{22})\delta \hfill \\ k_{n2}=\frac{\sqrt{3}}{3}(C_{11}+C_{22}+\sqrt{3}{C}_{16}+\sqrt{3}{C}_{26})\delta \hfill \\ k_{n3}=\frac{\sqrt{3}}{3}(C_{11}+C_{22}-\sqrt{3}{C}_{16}-\sqrt{3}{C}_{26})\delta \hfill \\ k_{s1}=\frac{2\sqrt{3}}{3}(3C_{66}-C_{22})\delta \hfill \\ k_{s2}=\frac{\sqrt{3}}{3}(C_{22}-3C_{12}+3\sqrt{3}{C}_{16}-\sqrt{3}{C}_{26})\delta \hfill \\ k_{s3}=\frac{\sqrt{3}}{3}(C_{22}-3C_{12}-3\sqrt{3}{C}_{16}+\sqrt{3}{C}_{26})\delta \hfill \end{array}$$

(1)where $k_{n1}$ and $k_{s1}$ are the normal and tangential spring constants between particle 0 and 1, and particle 0 and 4, respectively; $k_{n2}$ and $k_{s2}$ are the normal and tangential spring constants between particle 0 and 3, and particle 0 and 6, respectively; $k_{n3}$ and $k_{s3}$ are the normal and tangential spring constants between particle 0 and 2, and particle 0 and 5, respectively; δ is the diameter of particles in this study, which presents the thickness of a lamina. $C_{i,j}(i,j=1,2,6)$ are the elastic coefficients of the material stiffness matrix, which can be shown in equation (2):

$$C=\left(\begin{array}{ccc}{C}_{11}& C_{11}& C_{16}\\ C_{12}& C_{22}& C_{26}\\ C_{16}& C_{26}& C_{66}\end{array}\right)$$

(2)

The reduced stiffness matrix of the material can be obtained by the Young’s Modulus and Poisson’s ratios within the plane stress state when $C_{16}=C_{26}=0$ ,

$$\{\begin{array}{l}{C}_{11}=\frac{E_1}{1-{\nu }_{12}{\nu }_{21}}\hfill \\ C_{22}=\frac{E_2}{1-{\nu }_{12}{\nu }_{21}}\hfill \\ C_{12}=\frac{{\nu }_{12}{E}_2}{1-{\nu }_{12}{\nu }_{21}}\hfill \\ C_{66}=G_{12}\hfill \end{array}$$

(3)

For the 30∘ and 90∘ plies shown in Figure 2(b), the x and y axials in the original coordinate, E₁ and E₂, and ν₁₂ and ν₂₁ need to be exchanged with each other to obtain the bond normal and tangential stiffnesses. ³⁷

Modeling of $\pm 45°$ -ply lamina

Regular square packing is successfully applied to generate the model of $\pm 45°$ -ply lamina from our previous study and validated by the theoretical findings and FEM simulations. ³⁷ The configuration of the nine balls model and the model of $\pm 45°$ ply composite laminae can be found in Figure 4. $k_{n1}$ and $k_{s1}$ are the normal and tangential spring constants between particle 0 and 1, and particle 0 and 5, respectively; $k_{n2}$ and $k_{s2}$ are the normal and tangential spring constants between particle 0 and 2, between particle 0 and 4, between particle 0 and 6 and between particle 0 and 8, respectively; $k_{n3}$ and $k_{s3}$ are the normal and tangential spring constants between particle 0 and 3, and particle 0 and 7, respectively; δ is the particle diameter. For the anisotropic materials, k_n and k_s in equation (4) are calculated from the $C_{i,j}(i,j=1,2,6)$ , which are the elastic coefficients of the material stiffness matrix, which can be obtained from equation (3) under the plane stress state. The strength of bonds in $\pm 45°$ ply lamina can be calibrated from the strength predicted by the Tasi-Hill failure criterion, discussed in detail in the next section.

$$\{\begin{array}{l}{k}_{n1}=\frac{1}{3}(3C_{11}+C_{12}-4C_{66})\delta \hfill \\ k_{n2}=\frac{1}{3}(C_{12}+2C_{66})\delta \hfill \\ k_{n3}=\frac{1}{3}(C_{12}+3C_{22}-4C_{66})\delta \hfill \\ k_s=\frac{2}{3}(C_{66}-C_{12})\delta \hfill \end{array}$$

(4)

Figure 4.

Extended 2D square DEM model of $\pm 45°$ angle-ply composite lamina: (a) Configuration of square packing, (b) $\pm 45°$ ply composite lamina.

3D DEM models for cross-ply and QI CFRP composite laminates

Two experimental cases, one cross-ply laminates and one QI composite laminates with an open hole, were selected to validate and compare with the 3D DEM models in terms of cross-ply and angle-ply CFRP composite laminates. A carbon/epoxy composite material named Hexcel M21/T800 with a stacking sequence of ${[0_2^{°}/{90}_4^{°}]}_s$ was selected for the cross-ply laminates modeling via 3D DEM techniques. The specimens in the experiments ⁴⁰ were 180 mm long, 30 mm wide and 2.208 mm thick, and the number of cracks was counted within 30 mm of the specimen. Therefore, in order to reduce the computational cost, the length of 30 mm was adopted in the 3D DEM model and the width was also scaled down to one-sixth of the original one. In both cases, the thickness of each ply is equal to the diameter of particles, so that each ply constructed in the 3D DEM model can represent the realistic ply within the composite material systems. For the OHT QI composite laminates, a 16-ply ${[45°/0°/-45°/90°]}_{2s}$ IM7/8552 carbon/epoxy composite laminates test coupon ⁴¹ was selected, which is 190.5 mm long, 38.1 mm width and 2.642 mm thick, with a hole diameter of 6.35 mm. Considering the symmetry of geometry of the composite and computational effort, half of the coupon through the thickness is adopted in the modeling and symmetric boundary conditions were applied. Therefore, in this study, 33 052 particles are generated to create the 3D DEM model of the cross-ply laminate with 108 074 bonds, while the angle-ply laminate has 123 972 particles and 441 343 bonds.

When considering the damage progression of transverse cracking of composites, the strength distribution of the bonds with 90∘ plies should be taken into consideration. The accumulation of the transverse cracking results from their inherent defects such as microcracks, voids, areas of high fiber volume fraction, etc. causing the transverse layer to have a statistical distribution of strength along its length. ⁴² Thus, the random distribution of strength with a variation of 30% following a normal distribution law ⁴² was assigned in the 90∘ plies, while a variation of 10% was used in the strength distribution of bonds within 0∘ plies. A random generator was used in the DEM code PFC v5.30 to obtain the random numbers ranging within $(0.7-1.3)$ , which was assigned to the bonds of 90∘ plies by multiplying with the calibrated bond strength. Of note, in order to capture the fiber splitting in $\pm 45°$ plies in the QI composites, the strengths of fiber bonds and matrix bonds were assigned with the calibrated bond strength of 0∘ plies from experiments and the calibrated bond strength of $\pm 45°$ from the Tasi-Hill failure criterion. Then the cross-ply and QI CFRP laminates were subjected to the uniaxial tensile loading by applying a constant, low velocity on the red and black particles at both sides, shown as Figure 5. Figure 5 shows the 3D DEM models of ${[0_2^{°}/{90}_4^{°}]}_s$ M21/T800 CFRP cross-ply composite laminate and ${[45°/0°/-45°/90°]}_{2s}$ IM7/8552 CFRP composite laminate with an open hole. The particles are generated to create bonds, on which different mechanical properties can be applied. In both models, the particles, except those at the boundaries, are scaled down for a better view of the bonds at adjacent particles. The interface between different layers is created between two particles that share the same plane dimensions. It should be noted that in the QI composite laminate, the square pattern cannot match with the hexagonal one very well; an additional adjustment was made to use the rectangle pattern to approximate square pattern for $\pm 45°$ plies. Therefore, the number of interface bonds between $45°$ (or $-45°$ ) and 0∘ is twice as the bonds between $45°$ (or $-45°$ ) and 90∘, because of the inherent patterns of the different plies. It can be found in Figure 5 the cyan bonds of 0∘ plies in the longitudinal direction and the cyan bonds of 90∘ plies in transverse direction have the largest stiffness $1.2\times {10}^{15}\text{P}\text{a}/\text{m}$ calculated from the average strain energy, which are called fiber bonds. Purple bonds in 0∘ and 90∘ plies are called matrix bonds which have the normal stiffness of $4.39\times {10}^{13}\text{P}\text{a}/\text{m}$ . Red bonds between the neighbouring plies are regarded as interface bonds, which have the normal stiffness of $2.43\times {10}^{13}\text{P}\text{a}/\text{m}$ . It should be noted that the fiber direction cannot be presented in the current $\pm 45°$ plies by the normal stiffness of bonds; instead it is determined by assigning the fiber bond strength calibrated from 0∘ ply to the target bond ( $45°$ or $-45°$ ) to represent $45°$ or $-45°$ ply.

Figure 5.

3D DEM model for: (a) ${[0_2^{°}/{90}_4^{°}]}_s$ M21/T800 carbon/epoxy cross-ply composite laminates and (b) a 16-ply ${[45°/0°/-45°/90°]}_{2s}$ IM7/8552 carbon/epoxy laminate open-hole tensile test coupon.

The generation of each new crack can be captured when a bond breaks and various information can be recorded to evaluate the performance of the laminates during this loading process, such as crack density along the length (including the cracks in the 90∘ plies), the stress distribution of each ply, the stress-strain curve, the crack density-strain curve and the longitudinal Young’s modulus. The procedure is cycled until a given elongation is reached. The schematic diagram of the procedure is shown in Figure 6. In this model, the interfaces between different layers are modeled with the parallel bond model. The material properties of the unidirectional carbon fiber reinforced epoxy composite laminae can be found in Table 1.

Figure 6.

Schematic flowchart for the DEM modeling of damage progression.

Table 1.

Mechanical properties of UD carbon fiber epoxy composite laminae.

	M21/T800 40	IM7/8552 41
Longitudinal modulus E1 (GPa)	130.7	171
Transverse modulus E2 (GPa)	8.01	8.96
Poisson’s ratio ϵ12	0.35	0.32
Shear modulus G12 (GPa)	5.31	3.95
Transverse tensile strength YT (MPa)	40.5	98.6
Shear strength SL (MPa)	69.6	113

Failure of UD composite laminae

Tsai-Hill failure criterion

When simulating the failure process of the whole ply of composites, a failure criterion is usually needed in the numerical modeling with the framework of FEM. However, in the DEM-based modeling, due to the lack of rigorous formulation between the bond strength and macroscopic materials properties, the calibration of bond strength for any DEM model is inevitable. Bonds in 0∘ and 90∘ plies can be normally calibrated from the experimental ultimate tensile strength tests, while a failure criterion is needed for an angle-ply lamina when its experimental data are not available. Considering the simplicity of parameter identification of the criteria, Tsai-Hill failure criterion, ⁴³ which is based on the distorted energy theory and interaction failure theory, is adopted for this model. This criterion is a quadratic, interactive stress-based criterion that identifies failure but does not distinguish between different modes of failure. Although the Tsai-Hill failure criterion does not distinguish between the compressive and tensile strengths, which can result in an underestimate of the maximum loads that can be applied when compared to other failure theories like Tsai-Wu ⁴⁴ or Hashin ⁴⁵ theories, it can be still utilized in this study where only tensile loading condition is considered.

Considering a unidirectional composite material lamina is transversely isotropic, the Tsai-Hill failure criterion, which has been proved to be successful under a wide variety of circumstances, ⁴⁶ is given in equation (5) when the lamina is subjected to plane stress states.

$$\frac{{\sigma }_1^2}{X^2}-\frac{{\sigma }_1{\sigma }_2}{X^2}+\frac{{\sigma }_2^2}{Y^2}+\frac{{\tau }_{12}^2}{S^2}=1$$

(5)where X and Y are the tensile strength in x and y directions, respectively, S is the shear strength in x – y plane, and σ₁ and σ₂ are the principal stresses of the lamina in x and y directions.

For a unidirectional angle-ply (θ) lamina subjected to a uniaxial loading (σ_x), which can be found in Figure 7, the transformation equations for expressing stresses in a 1–2 coordinate system (principal stresses) in terms of stresses in an x–y coordinate system (applied stresses) is written as

$$\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\tau }_{12}\end{array}\right]={\left[T\right]}^{\hspace{0.17em}}\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]$$

(6)

Figure 7.

Failure strength prediction of a unidirectional lamina with Tsai-Hill failure criterion.

where

$$\left[T\right]=\left[\begin{array}{ccc}{\cos }^2\theta & {\sin }^2\theta & 2\cos \theta \sin \theta \\ {\sin }^2\theta & {\cos }^2\theta & -2\cos \theta \sin \theta \\ -\cos \theta \sin \theta & \cos \theta \sin \theta & {\cos }^2\theta -{\sin }^2\theta \end{array}\right]$$

(7)

When the lamina is under an applied uniaxial tension ( ${\sigma }_y={\tau }_{xy}=0$ ), then the principal stresses can be expressed by the applied stress as:

$$\{\begin{array}{l}{\sigma }_1={\sigma }_x{\cos }^2\theta \hfill \\ {\sigma }_2={\sigma }_x{\sin }^2\theta \hfill \\ {\tau }_{12}={\sigma }_x\sin \theta \cos \theta \hfill \end{array}$$

(8)combining equation (8) with equation (5), the Tsai-Hill failure criterion for the angle-ply lamina under uniaxial tension can be written as:

$$\frac{{\cos }^4\theta }{X^2}+\left[\frac{1}{S^2}-\frac{1}{X^2}\right]{\cos }^2\theta {\sin }^2\theta +\frac{{\sin }^4\theta }{Y^2}=\frac{1}{{\sigma }_x^2}$$

(9)

This criterion is used for the determination of failure strength of angle-ply composite laminae, which is then utilized to calibrate the strength of matrix bonds within the angle-ply laminae. Here in this study, only $45°$ plies are involved, so the calculated tensile strength of $45°$ ply σ_x from the Tsai-Hill failure criterion is adopted for the calibration of the matrix bond strength, while the fiber bond strength comes from the calibration of bond strength of 0∘ ply unidirectional composite laminae. For the sake of simplicity, it was assumed that the shear strength of bonds is equal to its normal strength. ³⁰ It is important to note that in the current study only tensile loading condition is considered to predict the damage initiation and progression in the composite laminates via the DEM model. The relationship between the ply angle and strength of angle-ply lamina of two different composite materials is plotted in Figure 7 with the Tsai-Hill failure criterion. It can be found that the strengths drop dramatically when the angle θ increases slowly before it reaches $10°$ , while it becomes stable when the angle is beyond $30°$ .

Definition and visualization of cracks

Figure 8 illustrates the schematic representation of cracks in a cross-ply/angle-ply laminate (up) and numerical representation of the cracks in 90∘ plies of a cross-ply composite laminate. Normally, a crack will be counted when the crack locates through the thickness, which can be shown in Figure 8(a1) and (a2). However, due to the randomness of the strength of bonds, a crack will be counted once the length of the consecutive or interrupted cracks is larger than half of the thickness of the weak plies, such as shown in Figure 8(b1) and (b2). Figure 8(b2) shows an example of the cracks in the 90∘ plies of a cross-ply laminate, as each ply is modeled separately within the composite, a red disk along the thickness direction represents a crack in this ply. The cracks can be initiated from the 90∘ plies ply adjacent to 0∘ plies at both side and propagate toward the middle plane, or from the 90∘ plies at the bottom (top) adjacent the 0∘ plies and propagate toward the top (bottom) plies in the direction perpendicular to the loading direction. Therefore, these two cases will be considered when counting the cracks during the loading.

Figure 8.

Schematic representation of cracks in a [ $0_m^{°}/{\theta }_n^{°}/0_m^{°}$ ] laminate (up) and numerical representation of 90∘ ply cracks in the ${[0_2^{°}/{90}_4^{°}]}_s$ composite laminates (down).

Case study of damage progression in composite laminates

Two case studies were conducted with this 3D DEM model by implementing the Tsai-Hill failure criterion for predicting the damage initiation and progression, including a cross-ply composite laminates and a QI composite laminates with a circular hole at the center.

Progressive failure in cross-ply laminates under tensile loading

In this study, an example of ${[0_2^{°}/{90}_4^{°}]}_s$ M21/T800 CFRP cross-ply composite laminates is modeled. Damages, such as transverse cracking and interfacial delamination in the laminates, are represented by the breakage of the bonds at the contacts between particles in the 3D DEM model. The bond at the interface between different layers plays a significant role by transferring the load from 0∘ ply to 90∘ ply, which is very sensitive to its shear stiffness and shear strength. In this study, the interlaminar shear strength is adopted from experimental data, ⁴⁷ while two different stiffnesses are considered to investigate the influences of the interfacial stiffness on the interlaminar delamination of composite laminates. Since the stiffness of the interface is difficult to be determined directly by experimental techniques and the use of a simple elastic brittle constitutive model for interlaminar interface, the bond stiffness of the interface is calculated based on the assumption that the bond stored strain energy (U) equals to the critical fracture energy release rate G_IIc. Here $U=G_{\mathit{IIc}}=\frac{1}{2}\tau {\delta }_s=\frac{{\tau }^2}{2\overline{k_s}}$ , then the bond shear stiffness is $\overline{k_s}$  =  $\frac{{\tau }^2}{2G_{\mathit{IIc}}}$ . Thus, the calculated shear stiffness of $1.4\times {10}^{12}\text{P}\text{a}/\text{m}$ and another larger value of $1.4\times {10}^{13}\text{P}\text{a}/\text{m}$ were adopted in this model to investigate the interlaminar delamination on the progressive damage within the composite laminates as earlier delamination interface cracks would be found in the DEM simulation with the larger interface shear stiffness according to the constitutive model of the parallel bond.

The images of sequential cracks in Figure 9(a1) to (a3) come from the specimen 1 in the experiment from Sebaey et al., ⁴⁰ where crack occurs since loading strain reaches 0.4%. However, in the other two specimens, the crack starts from 0.32%, which corresponds to the cracks obtained from the numerical simulations in Figure 9(a2) and (a3). Due to the experimental scattering data, the number of cracks in three specimens near the applied strain of 0.4% ranges from 2 to 5, which can be calculated from Figure 10. The cracks counted in the DEM modeling with the smaller stiffness of $1.4\times {10}^{12}\text{P}\text{a}/\text{m}$ (6 cracks) is 20% larger than those obtained from the experimental specimen 1 and 3, while the cracks in the DEM modeling with the higher stiffness of $1.4\times {10}^{13}\text{P}\text{a}/\text{m}$ (8 cracks) is 60% larger. It can be found from Figure 10 that more experimental data are closer to the numerical results obtained from the interface stiffness of $1.4\times {10}^{12}\text{P}\text{a}/\text{m}$ , which means this value can predict the progressive cracking better than the scaled-up one.

Figure 9.

Comparison of the results from numerical simulations and experiments ⁴⁰ regarding the sequent crack events of 90^∘ plies at different uniaxial strains: (a1,b1,c1) ${\epsilon }_x=0.3\%$ ; (a2,b2,c2) ${\epsilon }_x=0.4\%$ ; (a3,b3,c3) ${\epsilon }_x=0.5\%$ (a1,a2,a3 represent the experimental results; b1, b2, b3 represent the DEM simulation results with the smaller stiffness of $1.4\times {10}^{12}Pa/m$ ; c1, c2, c3 represent the DEM simulation results with the higher stiffness of $1.4\times {10}^{13}Pa/m$ ).

Figure 10.

Comparison of the results of crack density vs applied strain from experiments ⁴⁰ and DEM simulations, and stress-strain curve of simulations with and without interfacial delamination.

When the tensile stress is applied to cross-ply laminate along the fiber direction, cracks are formed in the 90∘ plies at a load level much lower than the ultimate strength of the laminate. ⁴⁸ In Figure 10 the small drop (bonds breakage) occurs at 0.22% in the DEM simulation with the interface bond stiffness of $1.4\times {10}^{13}\text{P}\text{a}/\text{m}$ and at 0.24% in the simulation with the stiffness of $1.4\times {10}^{12}\text{P}\text{a}/\text{m}$ , respectively. This slight difference exists due to the larger interfacial bond stiffness in the simulation with interfacial delamination, resulting in stress concentration in these weakest bonds, and these bonds only exist randomly inside the laminate so that the transverse cracks cannot be detected. Once the cracks are initialized from these small flaws in the 90∘ plies, they start to spread as tunneling cracks in the width direction (Figure 11(b)), and then span the entire 90∘ ply, but finally are arrested at the $0°/90°$ ply interfaces. ⁴⁹ However, the cracks in the DEM simulation with the larger stiffness (shown in Figure 11) suggest that the crack at the $0°/90°$ ply interfaces occurs around the strain of 0.3%. As the load increases further, more transverse cracks are initiated in the 90∘ plies and more interface cracks are formed, see Figure 11(c) to (g). It should be noted that the interesting phenomenon found in the DEM simulation with the smaller stiffness is that the interfacial cracks at $0°/90°$ ply were not observed during the whole loading process.

Figure 11.

Top view of sequent cracks of 90° plies from DEM simulation with interfacial delamination at different uniaxial strains in width direction (the pink disks represent the cracks at the $0°/90°$ ply interfaces). (a) Top view, (b) ε _x  = 0.3%, (c) ε _x  = 0.4%, (d) ε _x  = 0.5%, (e) ε _x  = 0.6%, (f) ε _x  = 0.7% and (g) (b) ε _x  = 0.8%.

The remaining stiffness of the laminates was calculated from the corresponding loading stress and strain at different loading strains and compared to the experimental and theoretical results in Sebaey et al., ⁴⁰ which are normalized in Figure 12. The same stiffness reduction rate can be observed between the experimental findings and numerical results before the applied strain of 0.3%, which corresponds to the same crack density at this strain found in Figure 10. However, a big difference is detected when the applied strain reaches 0.35% because the crack density obtained from the numerical simulation results with two different interface stiffnesses are much larger than that in the experiments. Moreover, the underestimated stiffness from the numerical simulations may be due to the fact that the experimental data were measured after the specimens were unloaded ⁵⁰ during which some micro-crack may be closed, resulting in a relatively higher remaining stiffness, while the stiffness was directly calculated by the stress and strain at the exact loading points in numerical simulations without a loading and reloading process. ³² In addition, the same trend for the curve of reduced stiffness can be found since the strain goes up to 0.5%, where the reduced stiffness reaches a plateau. Since it reaches the plateau, the slight difference between the numerical results with different interface stiffnesses suggests that the interfacial delamination has little influence on the remaining stiffness of the composite during the damage process, which can also be seen in Figure 10 where the two stress-strain curves overlap. This is mainly because the majority of the load is still held by the 0∘ plies. A small difference can be seen between the results from the DEM simulations and Ply Discount Method (PDM) since the remaining stiffness is calculated when the elastic properties (not include the Young Modulus in the fiber direction) of 90∘ plies are set to 0 in PDM. ⁴⁰ The damages in Figures 9(c3) and 11(g) suggest that the 90∘ plies can only hold little load due to the significant amount of cracks.

Figure 12.

Comparison of the results of reduced stiffness in longitudinal direction vs. applied strain from experiments, Ply Discount Method (PDM) ⁴⁰ and DEM models.

Progressive failure in the QI composite laminates with an open hole under tensile loading

In addition to the qualitative analysis for the QI ${[45°/0°/-45°/90°]}_{2s}$ IM7/8552 CFRP composite laminates with an open hole, a calibration of each ply is still necessary. Intralaminar bond strengths were calibrated based on the experimental data of longitudinal and transverse tensile strength ⁴¹ as well as the Tasi-Hill predicted strength of $\pm 45\%$ ply, while the interlaminar interface bond stiffness was calibrated from the experimental data.^41,51 Together with the maximum failure stress criterion, the damage initiation and propagation can be captured thanks to the natural characteristics of this DEM solution. It was reported in Hallett et al. ⁵² that the damage in OHT specimens can be split into four different stages according to its location, which is shown schematically in Figure 13, occurring in the following order:

Figure 13.

Damage initiation and propagation of ${[45°/0°/-45°/{90}^{°}]}_{2s}$ IM7/8552 CFRP composite laminates under tensile loading. ⁵²

1. Isolated crack/damage around the hole and free edge of the specimen;

2. Extensive damage around the hole (inner delamination regions);

3. Damage across the width of the specimen from the hole (outer delamination regions);

4. Final catastrophic failure.

Due to the nature of the randomness of bond strength in the composites, cracking occurs randomly in different plies at a relatively low strain, which is shown in Figures 14(a) and 15(a). Figure 14 shows the damage initiation and propagation in each ply of ${[45°/0^{°}/-45°/{90}^{°}]}_{2s}$ IM7/8552 CFRP composite laminates under tensile loading. It can still be observed that the first damage in the laminates occurs as matrix cracking in the 90∘ plies at the hole boundary perpendicular to the loading direction, which is followed by isolated cracks occurring in the off-axis plies within the inner delamination regions.

Figure 14.

Damage initiation and propagation of each ply within the ${[45°/0°/-45°/90°]}_{2s}$ IM7/8552 CFRP composite laminates under tensile loading at different strains: (a) ${\epsilon }_x=0.3\%$ , (b) ${\epsilon }_x=0.44\%$ , (c) ${\epsilon }_x=0.58\%$ and (d) ${\epsilon }_x=0.73\%$ . (The cracks are represented by deleting the bonds.).

Figure 15.

Damage initiation and propagation of interface of ${[45°/0°/-45°/90°]}_{2s}$ IM7/8552 CFRP composite laminates under tensile loading. (a) ε _x  = 0.3%, (b) ε _x  = 0.44%, (c) ε _x  = 0.58%, (d) ε _x  = 0.73%.

In the second stage of damage development, the isolated cracks and splits in different plies joined together with associated local delaminations at the inner delamination regions, as shown in Figure 15(b) and (c). The delamination was initiated at the hole edge nearest to the 90∘ plies and propagated along the thickness direction towards the middle plane through the cracks in the adjacent plies. However, when the delamination propagation stepping up and down reached the 0∘ plies, it would be obstructed by the 0∘ plies, thus changing the propagation to the length direction. It can be seen in Figures 14(b) and (c) and 15(b) and (c) that the cracks in 0∘ plies can be captured since the loading strain ε reaches 0.44%, and then propagate along fibers from the hole edge, resulting in the splits of fibers at the strain of 0.58%. Piles pulled apart from each other, and cracks in the different interfaces were found, where damage propagation in the form of delamination cracks passed through the neighboring interface. After the 90∘ plies lost the capability of carrying loads, the cracks propagated from the interface to $\pm 45°$ plies close to the middle plane, forming the splits in the $\pm 45°$ plies at the strain of 0.58%, see Figure 14(c).

In the third stage of damage development, the cracks propagated across the width from the inner delamination region to the outer delamination region until the free edge of the specimen. The first damage across the width of the specimen is in the 90∘ ply at the middle plane of the model, which is followed by $\pm 45°$ cracking close to the middle plane, as seen in Figure 14(c) and (d). When the $\pm 45°$ plies close to the middle plane lost their capability of carrying loads, the cracks initiated from the interface between 0∘ and $\pm 45°$ plies started to propagate into $\pm 45°$ plies near the surface, as shown in Figure 14(c). When the loading still increased, the cracks at the adjacent layers around the hole and free edge of the specimen started to propagate toward each other until they joined together, leading to full delamination across its width, as shown in Figure 15(d). Again, when the cracks started to propagate across the thickness, and reached the 0∘ plies, they changed the propagation direction toward the end of the specimen, which can be observed in the 0∘ plies.

The propagation process of the delamination is seen in Figure 15(c) and (d). In addition, when the stress concentration accumulated in the 0∘ plies, reaching the strength of the highly deformed bonds, and then the redistribution of stress would occur, thereby delaying the fiber fracture (see Figure 14(d)), as explained in Mollenhauer et al.. ⁵³ It should be noted that the delamination cracking is symmetric in the early loading stage while it becomes asymmetric when the strain is beyond 0.7%, which is probably due to the random distribution of bond strength at the interface, where the weakest bonds lead the progressive cracking path.

The damage zones in the third stage of damage are illustrated in Figure 16 with the comparison of results from experimental X-Ray techniques and DEM modeling. The delamination around the hole edge is clearly visible due to the high shear stress and Poisson’s ratio, and the same phenomenon is also found in the numerical results, where damage zones are represented by cracks at different loading strains (Figure 16(b) and (c)). Figure 16(b) suggests that the symmetric damage zones represented by the superimposed cracks start to emerge when the loading strain reaches 0.58%. As the strain increases until 0.73%, the triangular damage zone becomes more obvious, as shown in Figure 16(c), in which the yellow lines follow the edge of the damage regions.

Figure 16.

Damage zones associated with ${[45°/0°/-45°/90°]}_{2s}$ IM7/8552 CFRP composite specimens: (a) X-Ray experimental result, ⁴¹ (b) damage zone represented by cracks at strain of 0.58%, (c) damage zone represented by cracks at strain of 0.73%.

The final catastrophic failure occurred in the form of fiber fracture when loading stress reached 378 MPa in the experiments and 340 MPa in the numerical simulation. Figure 17 shows the comparison of final surface-ply ( $45°$ ply) cracks of the experiments and the numerical simulation in the ${[45°/0°/-45°/90°]}_{2s}$ IM7/8552 CFRP composite laminates under tensile loading. After the specimen experienced the three stages of damage development, it reached the failure strength of the fiber within 0∘ plies. It can be found in the comparison that the DEM simulation result shows a good agreement with the two experimental findings, together with the fiber fracture shown in the Figure 14(d) close to the middle plane, suggesting that the DEM model is capable of capturing the final failure of the whole composite laminates.

Figure 17.

Comparison of final surface-ply cracks in the ${[45°/0°/-45°/90°]}_{2s}$ IM7/8552 CFRP composite laminates experimental OHT coupons (σ_x = 378 MPa) ⁴¹ and numerical simulation result (σ_x = 360 MPa) under tensile loading: (a) OHT1 (b) OHT2 (c) DEM model.

Conclusions

In this study, the progressive failure of carbon fiber-reinforced composite laminates under uniaxial tension was numerically investigated by means of 3D DEM modeling. Extended 2D hexagonal and square packing strategies were adopted for the construction of 0∘ and 90∘ plies, and 45∘ ply, respectively. The parallel bond model was utilized to capture the constitutive behaviors of the bond between neighboring particles. The relationship between the bond stiffness and material properties was established based on the average strain energy method. The random distribution of bond strengths calibrated from experiments with a variation of 30% and 10% following a normal distribution law are assigned to the bonds in 0∘ and 90∘ plies to capture random cracks, respectively. The bond strength of 45% ply in the composite was calibrated from the failure strength prediction by the Tsai-Hill failure criterion. The maximum stress failure criterion was applied to predict the bond breakage of the DEM models.

Quantitative analysis was conducted for progressive damage prediction of the cross-ply M21/T800 CFRP with a stacking sequence of ${[0_2^{°}/{90}_4^{°}]}_s$ . Two different interface shear stiffnesses were utilized for the prediction of cracking in the composite. It was found that the interface stiffness has little effect on the stress-strain curve, which is mainly because the 0∘ plies still hold the load after cracking in the 90∘ plies and interface occurs. However, it has significant influences on the crack density in the 90∘ plies, since the larger interface stiffness can transfer the load faster than the one calculated from fracture energy between neighboring plies, resulting in a higher crack density. The value of interface bond stiffness calculated from the fracture energy is recommended due to its theoretical foundation, and better prediction can be achieved with a more sophisticated bond model.

The qualitative analysis of failure prediction of the ${[45°/0°/-45°/90°]}_{2s}$ QI IM7/8552 carbon/epoxy composite laminates with an open hole subjected to tensile loading demonstrates that the 3D DEM model is capable of modeling all the four stages of damage development as observed in experiments, including isolated cracks, inner delamination cracks, outer delamination cracks and final failure. Inner damage zones and final failure of the surface ply can be found in a good quantitative agreement with the experimental CT findings.

Besides, it is important to note that this is the preliminary attempt to investigate the damage initiation and propagation using 3D DEM solution. Currently, we are working on the quantitative analysis of the OHT problem with our experimental data. With the implementation of a more sophisticated cohesive model or seven-bond interface model for the interlaminar interface, the developed 3D model could be applied to predict the damage progression and failure strength of available composite laminates under more complex multiaxial loading conditions. From the perspective of computational cost, the simulations of cross-ply composite laminates under tensile loading cost around 9 h to reach its prescribed loading strain, while the simulation of the QI composite laminates with an open hole costs 40 h to reach its final failure. Of note, the automatic timestep mode should be utilized since it can adjust automatically and accordingly when bonds break, based on the current stiffnesses and masses of all objects in the system. In the future, GPU-based DEM code should be developed for the automatic generation of composite lay-up process and saving the computational cost, when more detailed damages in composite laminates can be captured with a larger DEM model.