# A Delamination Propagation Model for Fiber Reinforced Laminated Composite Materials

Authors: David Aveiga; Marcelo L. Ribeiro
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The employment of composite materials in the aerospace industry has been gradually considered due to the fundamental lightweight and strength characteristics that this type of materials has. The science material and technological progress reached matched perfectly with the requirements for high-performance materials in aircraft and aerospace structures; thus, the development of primary structure elements applying composite materials became something very convenient. It is extremely important to pay attention to the failure modes that influence composite materials performances, since these failures lead to a loss of stiffness and strength of the laminate. Delamination is a failure mode present in most of the damaged structures and can be ruinous, considering that the evolution of interlaminar defects can carry the structure to a total failure followed by its collapse. The present work aims at the development of a delamination propagation model to estimate a progressive interlaminar delamination failure in laminated composite materials and to allow the prediction of material’s degradation due to delamination phenomenon. Experimental data, available at literature, was considered to determine some model parameters, like the strain energy release rate, using GFRPs laminated composites. This new delamination propagation model was implemented as subroutines in FORTRAN language (UMAT-User Material Subroutine) with formulations based on the Fracture Mechanics and Continuum Damage Mechanics. Finally, the UMAT subroutine was complemented with an intralaminar model and compiled beside the commercial Finite Element (FE) software ABAQUS™.


---

## Body

## 1. Introduction

Composite materials afford the unique possibility of designing the material, the manufacturing process, and the structure in one unified and concurrent procedure; having a large number of degrees of freedom available enables simultaneous material optimization for several given constraints. The applications of composites tend to became highly desirable in primary and secondary aircraft structures, for example, wing stringers, floor beams, fuselage or wing skins surfaces, and the empennage.Along with the high strength characteristic, lightweight is a crucial attribute to consider composite materials as a sophisticated and high-performance material. Especially for aircraft operators, weight saving may be a significant aspect of decreasing the operational expenses since this action is related to combustible consumption. All around the world, aircraft manufactures took advantage and started applying composites to their aircraft structures as is noticed in modern airplanes like Boeing’s B787 or Airbus’s A380.One of the peculiarities of composite materials is the complexity of failure behavior due to the presence of three different material phases (fibers, matrix, and their interface) that compose this material and directly affect damage mechanisms. With safety being one of the mainstays of aviation, all this versatility of composite materials is limited, since failure within primary structures may produce catastrophic accidents. Engineers and structure designers must use large safety factors in order to guarantee an invulnerable performance of a composite part during its life, which leads to an overweight penalty and compromises the characteristic performance-weight ratio of laminated composite materials. Understanding and describing the failure modes that affect composites materials became substantial in order to develop more confident and efficient composite structures.Material scientists have proposed several approaches to characterize damage within composite but the discussion is still open. For a better understanding, the subject has been divided into two subsections: intralaminar failure regarding damage within the lamina like matrix microcrack or fiber breakage and interlaminar failure like delamination as a product of intralaminar damage propagation.Ribeiro, Tita, and Vandepitte [1] proposed an intralaminar model that estimates the material degradation due to matrix and/or fiber failure for tensile or compression loads. The present work aimed at the complementation of Ribeiro’s material model by describing the delamination propagation as a means of interlaminar fracture failure. Values for strain energy released rates (SERR) were obtained from fracture tests assisted by ASTM standards for Mode I, Mode II, and Mixed-Mode interlaminar fracture of laminated composite materials.With basis in Fracture Mechanics (FM), a Cohesive Law was completely formulated with the employment of a polynomial cubic function. Using Continuum Damage Mechanics (CDM), damage variables were computed to represent material degradation. The complete model was finally implemented as a subroutine (User Material-UMAT) in FORTRAN language and compiled beside the commercial FE software ABAQUS.Further sections will introduce the damage model by describing the theories considered, the formulations, and the numerical implementation. Finally, a simulation of a 3-point bending tests is performed and compared with experimental results in order to demonstrate the potentiality of the model proposed.


## 2. Basis for the Interlaminar Model

### 2.1. Cohesive Zone Model Approach

The degradation of the material or the damage evolution is controlled by a damage variabledas known from CDM; the interpretation of damage starts from the formulation of a Cohesive Law or Traction-Separation Law (TSL). A TSL, in Finite Element Analysis, represents a contact behavior between two surfaces and is strongly related to energy dissipation due to crack propagation (surface separation). To completely formulate a TSL, three specific cohesive parameters are needed: critical cohesive strengths, fracture toughness, and initial stiffness [2].Within Finite Element Analysis, the initial stiffness is considered as a suitable fixed value of 106 N/mm3, since this value helps with numerical convergence, and the diversity of other values, experimentally calculated, does not imply a significant difference of the results [3].The fracture within composite materials has a plastic-elastic behavior; thus, the fracture toughness is measured in terms of strain energy released rate, denoted byG with unit of J/m2. These rates of energy dissipated are calculated for each fracture mode (Modes I, II, and III) and for Mixed-Modes by performing experimental or numerical tests specific for each mode.The critical cohesive strengths or nominal stresses are calculated from values of SERR with the employment of the Energy-Balance approach developed by Griffith [4] and adapted by Irwin [5] and Orowan [6] for more ductile materials like composites: (1)σC=EGCπawhereσc is the nominal stress,GC is the SERR,a is the delamination length, andEis the laminate propertyE33 since the crack is not penetrating the fibers but the matrix.With these three parameters defined, a TSL can be formulated. A TSL is governed by a function that describes the dependence of a tractionT on a separation distance δ, as a function Tδ. By this, there are some statements to consider:(i)
A parameter represents the maximum traction sustainable by an element. It is denoted asTC.(ii)
Elements totally fail at a maximum opening. This is represented asδ0.(iii)
The area under the TSL curve must be equal to the critical SERR at each mode.It is widely accepted that polymeric based composites have a ductile behavior in the presence of plastic yielding. Specifically when speaking about fracture, it is usually observed that the crack grows within the matrix with a progressive advance and with local crack tip fields affected by fiber orientation. This behavior often requires a nonlinear function of the TSL in order to well represent the degradation. For present case study, a polynomial cubic function (see Figure1) is proposed.Figure 1
TSL with a polynomial cubic function.


### 2.2. Model Formulation

After reviewing a few TSL shape comparisons for composite materials in the literature [7, 8] and after testing different functions, a cubic polynomial function was selected. The TSL is governed by the following function: (2)Tδ=φTCδδ01-δδ02whereTδ is any traction or cohesive strength at the curve depending on a separation δ,φ is an adjustment coefficient,TC represents the critical cohesive strength for the mode selected, δ is the separation, and δ0 is the ultimate separation.Regardless of its shape, the area under the traction-separation curve must be equal to the corresponding critical SERR; this characteristic is earned from Griffith’s Theory. This leads us to find a value forδ0 by satisfying the following condition: (3)∫0δ0Tδdδ=GCThus,(4)δ0=4GCφTCThe critical value of the separation, where the maximum cohesive strength is located, can be found by setting the derivative of the function to zero:(5)dTδdδ=0Thus,(6)δC=δ03If we considerT=TC, then δ=δC. Solving for (2) will give us a value for φ: (7)φ≈2.6=135By following this procedure, a polynomial cubic TSL can be formulated for cases involving laminated composites with polymeric matrixes.The damage provoked by delamination will be controlled by damage variables related to energy dissipation. Within CDM, a damage variable is defined as the ratio between a damaged state and an undamaged state associated with local discretization. For our case study, the area delimited by the traction-separation function and a line drawn from the origin to any degraded point at the curve is interpreted as a damaged state (GD), while the whole area under the traction-separation curve is interpreted as the undamaged state (GC) (see Figure 2) [9].Figure 2
Damaged and undamaged areas.Thus,(8)d=AreaDamagedAreaUndamaged=GDGCThe formulation of TSLs and the determination of the damage variables must be performed for Mode I, Mode II, and Mixed-Mode interlaminar failure in order to totally characterize the delamination propagation. The energy dissipation at Mode III fracture has been demonstrated to be insignificant when compared to the total energy released in a multiaxial load case [10]. This evidence is based on the fact that the shear in any direction, within the crack plane between the fracture surfaces, is the same under subsequent Mode II or Mode III loading. This study will consider the energies released rates GII and GIII as equivalent.


## 2.1. Cohesive Zone Model Approach

The degradation of the material or the damage evolution is controlled by a damage variabledas known from CDM; the interpretation of damage starts from the formulation of a Cohesive Law or Traction-Separation Law (TSL). A TSL, in Finite Element Analysis, represents a contact behavior between two surfaces and is strongly related to energy dissipation due to crack propagation (surface separation). To completely formulate a TSL, three specific cohesive parameters are needed: critical cohesive strengths, fracture toughness, and initial stiffness [2].Within Finite Element Analysis, the initial stiffness is considered as a suitable fixed value of 106 N/mm3, since this value helps with numerical convergence, and the diversity of other values, experimentally calculated, does not imply a significant difference of the results [3].The fracture within composite materials has a plastic-elastic behavior; thus, the fracture toughness is measured in terms of strain energy released rate, denoted byG with unit of J/m2. These rates of energy dissipated are calculated for each fracture mode (Modes I, II, and III) and for Mixed-Modes by performing experimental or numerical tests specific for each mode.The critical cohesive strengths or nominal stresses are calculated from values of SERR with the employment of the Energy-Balance approach developed by Griffith [4] and adapted by Irwin [5] and Orowan [6] for more ductile materials like composites: (1)σC=EGCπawhereσc is the nominal stress,GC is the SERR,a is the delamination length, andEis the laminate propertyE33 since the crack is not penetrating the fibers but the matrix.With these three parameters defined, a TSL can be formulated. A TSL is governed by a function that describes the dependence of a tractionT on a separation distance δ, as a function Tδ. By this, there are some statements to consider:(i)
A parameter represents the maximum traction sustainable by an element. It is denoted asTC.(ii)
Elements totally fail at a maximum opening. This is represented asδ0.(iii)
The area under the TSL curve must be equal to the critical SERR at each mode.It is widely accepted that polymeric based composites have a ductile behavior in the presence of plastic yielding. Specifically when speaking about fracture, it is usually observed that the crack grows within the matrix with a progressive advance and with local crack tip fields affected by fiber orientation. This behavior often requires a nonlinear function of the TSL in order to well represent the degradation. For present case study, a polynomial cubic function (see Figure1) is proposed.Figure 1
TSL with a polynomial cubic function.


## 2.2. Model Formulation

After reviewing a few TSL shape comparisons for composite materials in the literature [7, 8] and after testing different functions, a cubic polynomial function was selected. The TSL is governed by the following function: (2)Tδ=φTCδδ01-δδ02whereTδ is any traction or cohesive strength at the curve depending on a separation δ,φ is an adjustment coefficient,TC represents the critical cohesive strength for the mode selected, δ is the separation, and δ0 is the ultimate separation.Regardless of its shape, the area under the traction-separation curve must be equal to the corresponding critical SERR; this characteristic is earned from Griffith’s Theory. This leads us to find a value forδ0 by satisfying the following condition: (3)∫0δ0Tδdδ=GCThus,(4)δ0=4GCφTCThe critical value of the separation, where the maximum cohesive strength is located, can be found by setting the derivative of the function to zero:(5)dTδdδ=0Thus,(6)δC=δ03If we considerT=TC, then δ=δC. Solving for (2) will give us a value for φ: (7)φ≈2.6=135By following this procedure, a polynomial cubic TSL can be formulated for cases involving laminated composites with polymeric matrixes.The damage provoked by delamination will be controlled by damage variables related to energy dissipation. Within CDM, a damage variable is defined as the ratio between a damaged state and an undamaged state associated with local discretization. For our case study, the area delimited by the traction-separation function and a line drawn from the origin to any degraded point at the curve is interpreted as a damaged state (GD), while the whole area under the traction-separation curve is interpreted as the undamaged state (GC) (see Figure 2) [9].Figure 2
Damaged and undamaged areas.Thus,(8)d=AreaDamagedAreaUndamaged=GDGCThe formulation of TSLs and the determination of the damage variables must be performed for Mode I, Mode II, and Mixed-Mode interlaminar failure in order to totally characterize the delamination propagation. The energy dissipation at Mode III fracture has been demonstrated to be insignificant when compared to the total energy released in a multiaxial load case [10]. This evidence is based on the fact that the shear in any direction, within the crack plane between the fracture surfaces, is the same under subsequent Mode II or Mode III loading. This study will consider the energies released rates GII and GIII as equivalent.


## 3. Intralaminar Model

Ribeiro, Tita, and Vandepitte [1] proposed an intralaminar failure model regarding composite lamina under plane stress and considering uniform damage through the thickness. The model has three damage variables, each of them related to fiber failure, matrix failure, and shear stress failure under tensile or compression loads, respectively.Fiber failure is considered as linear elastic one with brittle fracture. At this direction, the model considers that matrix failure does not imply a degradation ofE11 property since the whole load is supported by the fibers. For this type of failure, an internal damage variabled1 was assigned.Matrix failure is exclusively controlled by transverse stress (σ22) and shear stress (τ12). At this direction, nonlinear behavior is observed due to inelastic strains [11]. A corresponding damage variable was assigned for each stress:d2 for transverse loading andd6 for shear loading. The hypotheses of effective stress were employed to link damage variables to the stresses: (9)σ˙11σ˙22σ˙12=11-d100011-d200011-d6σ11σ22τ12Table1 resumes the complete damage model and a full description can be found in [1]. The model has shown good potentiality by predicting progressive failure when loading has predominating in-plane stresses. Some application can be found in literature like a 3-point bending simulation of flat filament wound laminates [12], optimization of composite stacking sequence [13], hydrostatic external pressure of composite tubes [14], and radial compression [15].Table 1
Ribeiro's intralaminar failure model.
Failure Criteria Mode of Failure Degradation Law σ 11 X T ≤ 1 Fiber tensile E 11 = 0 σ 11 X C 0 ≤ 1 Fiber compression E 11 = X C 0 ε 11 1 - f ε 11 + f ε 11 E 1 1 0 f ≥ 0 Matrix tensile d 2 = A θ Y 2 + B θ f ≥ 0 Matrix compression E 22 = σ 2 2 y ε 22 1 - f ε 22 + f ε 22 E 2 2 0 f ≥ 0 Shear d 6 = C θ Y 6 + D θHeref is defined in (10) andf(ε11) is obtained from the fitting of stress-strain data of specimens under compressive loads. (10)f=σ222+τ122--S12y+2S12y1+σ22/σ2203Despite the good predictions achieved by employing this model, a reinforcement is required in order to improve the model predictions for cases when out-of-plane shear stresses are predominant.


## 4. Interlaminar Damage Model

A new damage variable is defined asd3 and is totally related to delamination phenomenon. It will affect the out-of-plane stress state of the material stiffness matrix. The value of this variable depends on the mode of fracture being considered. Table 2 presents the critical values for Mode I [16], Mode II [17], and Mixed-Mode [18] of carbon-epoxy material. For Mixed-Mode three different mode mixtures (GII/GC) were considered.Table 2
Critical values of SERR for carbon-epoxy.
Mode I Mode II Mixed-Mode G I I / G C = 20% G I I / G C = 50% G I I / G C = 80 % Carbon-Epoxy (kJ/m2) 1.439 2.11 0.28 0.43 0.62Using these values, five TSLs (Table3) can be generated by following the procedure described in Section 2.1. The maximum traction or critical nominal stress for each fracture mode is estimated with (1). From the plots, the damage variables are calculated by employing (8).Table 3
Nominal stresses and critical damage variable value.
Mode I Mode II Mixed-Mode G I I / G C = 20% G I I / G C = 50% G I I / G C = 80 % Nominal Stress (MPa) 35.81 55.98 6.26 7.76 9.32  Damage Variable (d 3 ) 0.87 0.91 0.89 0.89 0.89The interlaminar failure within carbon fiber composites has a brittle behavior; thus, the damage variables have almost the same critical value.A quadratic least square regression was performed using the values of damage variables for each fracture mode considered; this generates a single path for the variable evolution during simulations. The criterion also depends on the dominant direction of the loads. The next section gives details about model implementation via User Material Subroutine.


## 5. UMAT Implementation

The full model implementation considers two tridimensional user materials, one for intralaminar damage and another for delamination. The intralaminar damage model will degrade in-plane elastic properties only; on the other hand, the interlaminar model will degrade stiffness parameters.The stiffness matrix for intralaminar damage model is considered orthotropic:(11)D=1-υ23υ32E22E33Δ1-d2υ21+υ23υ31E22E33Δ1-d2υ31+υ21υ32E22E33Δ1-d2000υ21+υ23υ31E22E33Δ1-d21-υ13υ31E11E33Δ1-d1υ32+υ31υ32E11E33Δ1-d1000υ31+υ21υ32E22E33Δ1-d2υ32+υ31υ32E11E33Δ1-d11-υ21υ12E22E11Δ3000000G121-d6000000G13000000G23where(12)Δ=1-υ12υ21-υ23υ32-2υ12υ13υ21E11E22E331-d11-d2And(13)Δ3=1-υ12υ21-υ23υ32-2υ12υ13υ21E11E22E33The stiffness matrix for interlaminar damage model is considered isotropic:(14)Cijkl=λ∗δijδkl+μ∗δikδjl+υδilδjkwhere(15)μ∗=E∗21+υAnd(16)λ∗=E∗υ1+υ1-2υwhereυ is Poisson’s ratio and E∗=E(1-d3) withE being Young’s modulus.The damage variables introduced in the whole model are the maximum calculated values along the load history analyses in order to avoid material self-healing.


## 6. Experimental Testing and Numerical Simulation

### 6.1. Experimental Test

A 3-point bending test was carried out employing five specimens manufactured using an infusion process of preimpregnated carbon fibers. The specimen geometry and stacking sequence are given in Table4, and Table 5 shows the mechanical properties for the material considered.Table 4
Specimen geometry and stacking sequence.
Length Width Thickness Stacking Sequence (mm) (mm) (mm) Carbon/epoxy 70 10 3 0 ° 8Table 5
Material mechanical properties.
E 11 E 22 ν 12 G 12 G 23 (GPa) (GPa) (GPa) (GPa) Carbon/epoxy 140 11 0.3 4.5 2.2These properties were calculated at the laboratory of the Group of Aeronautical Structures (GEA from Portuguese) of the University of São Paulo localized in São Carlos City, excepting the value for out-of-plane shear modulus at the 2-3 planes (G23), which correspond to 20% of E2 [19].The test was performed employing an INSTRON universal testing machine adapted for 3-point bending with dimensions illustrated in Figure3. The test was displacement-controlled with a displacement rate of 0.5 mm/min.Figure 3
3-point bending test.


### 6.2. Numerical Simulation

Commercial FE software ABAQUS was employed to simulate the 3-point bending test following the same specifications detailed in Section6.1. A 3D solid deformable part was modeled (Figure 4) with the same geometry and material properties. The mesh had 53200 quadratic hexahedral elements of type C3D20, which are general purpose quadratic brick elements with three integration points, excellent for linear elastic calculations.Figure 4
Numerical simulation.The laminate stacking sequence was modeled by creating partitions for each ply and for each interface between the plies (see Figure5). Two different materials were defined for the lamina and for the interface between laminae, respectively. We understand these interface plies as representing an epoxy region between fiber laminates [20]; therefore, they are modeled as an isotropic material.Figure 5
Laminate design.The failure criteria were implemented in the UMAT subroutines using FORTRAN language in order to capture both intra- and interlaminar failure while running the simulation.The model presents almost no convergence problem and a promissory result discussed in the further section.


## 6.1. Experimental Test

A 3-point bending test was carried out employing five specimens manufactured using an infusion process of preimpregnated carbon fibers. The specimen geometry and stacking sequence are given in Table4, and Table 5 shows the mechanical properties for the material considered.Table 4
Specimen geometry and stacking sequence.
Length Width Thickness Stacking Sequence (mm) (mm) (mm) Carbon/epoxy 70 10 3 0 ° 8Table 5
Material mechanical properties.
E 11 E 22 ν 12 G 12 G 23 (GPa) (GPa) (GPa) (GPa) Carbon/epoxy 140 11 0.3 4.5 2.2These properties were calculated at the laboratory of the Group of Aeronautical Structures (GEA from Portuguese) of the University of São Paulo localized in São Carlos City, excepting the value for out-of-plane shear modulus at the 2-3 planes (G23), which correspond to 20% of E2 [19].The test was performed employing an INSTRON universal testing machine adapted for 3-point bending with dimensions illustrated in Figure3. The test was displacement-controlled with a displacement rate of 0.5 mm/min.Figure 3
3-point bending test.


## 6.2. Numerical Simulation

Commercial FE software ABAQUS was employed to simulate the 3-point bending test following the same specifications detailed in Section6.1. A 3D solid deformable part was modeled (Figure 4) with the same geometry and material properties. The mesh had 53200 quadratic hexahedral elements of type C3D20, which are general purpose quadratic brick elements with three integration points, excellent for linear elastic calculations.Figure 4
Numerical simulation.The laminate stacking sequence was modeled by creating partitions for each ply and for each interface between the plies (see Figure5). Two different materials were defined for the lamina and for the interface between laminae, respectively. We understand these interface plies as representing an epoxy region between fiber laminates [20]; therefore, they are modeled as an isotropic material.Figure 5
Laminate design.The failure criteria were implemented in the UMAT subroutines using FORTRAN language in order to capture both intra- and interlaminar failure while running the simulation.The model presents almost no convergence problem and a promissory result discussed in the further section.


## 7. Results and Discussion

The numerical simulation, working together with the UMAT, was able to reproduce the maximum load and, respectively, load displacement. Table6 resumes the maximum load and the respective displacement of each specimen tested and simulation. Figure 6 shows the load-displacement curve obtained.Table 6
Maximum load and respective displacement, experimental and numerical.
01 02 03 04 05 mean num error Maximum Load (kN) 1.52 1.57 1.29 1.41 1.34 1,43 1.64 13%  Load Displacement (mm) 2.41 2.51 2.03 2.36 2.26 2.31 2.26 2.4%Figure 6
Load-displacement curve.As seen in Figure6, the model was able to capture degradation due to interlaminar failure. The evolution of the delamination damage variable (d3), along the simulation, was captured for the nearest interface ply to the loading roller and is shown in Figure 7.Figure 7
Evolution of the delamination damage variable (d3).The results of maximum load and respective load displacement are acceptable. At the end of the first portion of the curves, before the maximum load, all the experimental specimens present a yield zone affected by matrix crushing (very common during 3-point bending test), which the numerical model could not represent. This matrix effect needs to be studied and characterized in order to be incorporated within the damage model. Another deviation that increases the discrepancy is the fact of considering a perfect specimen geometry during the simulation, where that the manufacture process produces specimens with variable geometry. The presence of diverse microdefects is not considered within the modeling.


## 8. Conclusions

A new delamination propagation model was introduced and showed good prediction of interlaminar failure. It works well in complement with Ribeiro’s intralaminar model to capture degradation of composite materials. The implementation of the model, using FORTRAN UMAT subroutine, is simple and demands low computational cost compared with other methodologies. The accuracy of the result can be considered acceptable.The model was able to capture delamination propagation and material degradation due to interlaminar failure. We can conclude that the proposed model has a strong potentiality to simulate and prevent delamination effects within composite structures.Some discrepancies were observed in the results due to the incapability of the model to capture matrix crushing which is present in bending loadings and the disregard of some material details. Additional studies are needed to account for these effects and may be proposed in further investigations.



---

*Source: 1861268-2018-06-19.xml*