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™.
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 variable d as 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 by G 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:
where σc is the nominal stress, GC is the SERR, a is the delamination length, and E is the laminate property E33 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 traction T on a separation distance δ, as a function . By this, there are some statements to consider:(i)A parameter represents the maximum traction sustainable by an element. It is denoted as .(ii)Elements totally fail at a maximum opening. This is represented as .(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 Figure 1) is proposed.
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:
where 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 by satisfying the following condition:
Thus,
The critical value of the separation, where the maximum cohesive strength is located, can be found by setting the derivative of the function to zero:
Thus,
If we consider , then . Solving for (2) will give us a value for :
By 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].
Thus,
The 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 and 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 of E11 property since the whole load is supported by the fibers. For this type of failure, an internal damage variable d1 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 and d6 for shear loading. The hypotheses of effective stress were employed to link damage variables to the stresses:
Table 1 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].
Here f is defined in (10) and f(ε11) is obtained from the fitting of stress-strain data of specimens under compressive loads.
Despite 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 as d3 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 () were considered.
Using these values, five TSLs (Table 3) 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).
The 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:
where
And
The stiffness matrix for interlaminar damage model is considered isotropic:
where
And
where is Poisson’s ratio and =E(1-d3) with E 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 Table 4, and Table 5 shows the mechanical properties for the material considered.
These 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 Figure 3. The test was displacement-controlled with a displacement rate of 0.5 mm/min.
6.2. Numerical Simulation
Commercial FE software ABAQUS was employed to simulate the 3-point bending test following the same specifications detailed in Section 6.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.
The laminate stacking sequence was modeled by creating partitions for each ply and for each interface between the plies (see Figure 5). 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.
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. Table 6 resumes the maximum load and the respective displacement of each specimen tested and simulation. Figure 6 shows the load-displacement curve obtained.
As seen in Figure 6, 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.
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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Disclosure
This manuscript is based on the partial conclusions of the research thesis of the author David Aveiga.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The support of Professor Marcelo L. Ribeiro, who provided insight and expertise that greatly improved the outcome of this study, is acknowledged. The authors acknowledge the São Paulo Research Foundation [FAPESP 2015/13844-8] and CAPES Foundation [nos. 1565501, PROEX-0487] for the financial support.