Cohesive Zone Modeling

NASTRAN COHESIVE ZONE MODELING OF COMPOSITE BONDED JOINTS Allan Hsu, Barton McPheeters NEi Software 5555 Garden Grove Boulevard Westminster, CA 92683, USA

ABSTRACT A computationally efficient method for modeling delamination initiation and propagation in composite bonded joints is presented to aid in aircraft structural analysis. The modeling strategy incorporates the following state of the art finite element technologies to accurately capture the complex damage mechanisms in bonded composite joints: cohesive zone modeling (CZM), three-dimensional layered solid elements, and micromechanics-based damage degradation. Simulation results are presented for a Double Cantilevered Beam (DCB) Mode I Fracture Test and a Mixed Mode Skin Stiffener Tension Test, and compared to experimental data from technical publications.

1. INTRODUCTION Analytical assessment of composite bonded joints is a recurring issue in aircraft design today. Various techniques have been applied to this problem with mixed results. Some techniques involve a two-step approach where a conventional finite element stress analysis predicts delamination initiation, and the initial debond results are passed to a fracture mechanics approach for delamination propagation predictions. Other techniques are not able to fully capture all of the critical components of complex damage mechanisms: initial damage resulting from through thickness matrix cracking and branching, matrix cracks driving delaminations to occur, stiffness degradations in fiber and/or matrix directions, and complete delamination across the adhesive bondline. This particular one-step modeling strategy is designed to capture the true failure mechanism of the joint, and incorporates the following state of the art finite element technologies to accurately capture the complex damage mechanisms in bonded composite joints: cohesive zone modeling (CZM), three-dimensional layered solid elements, and micromechanics-based damage degradation. As the micromechanics-based failure model, MultiContinuum Theory (MCT) is utilized as part of the modeling strategy to extract separate failure predictions for the fibers and epoxy matrices throughout the composite joints. Using a Double Cantilevered Beam (DCB) Mode I Fracture Test and a Mixed Mode Skin Stiffener Tension Test, the paper compares critical delamination propagation failure loads predicted by NASTRAN with experimental data from previously published technical publications.

1

http://www.nenastran.com/Composites

1. FINITE ELEMENT METHODS Cohesive Zone Modeling (CZM) is a generally applicable nonlinear methodology for simulating crack initiation and crack growth in pure mode and mixed mode fractures. The incremental growth of a crack due to environmental loading can be monitored from analysis, and the damage tolerance of the composite structure can be assessed. Many different delamination scenarios can be analyzed with CZM: delamination in the adhesive/adherend interface of a bonded joint or delamination in the skin/core interface of a large scale honeycomb sandwich structure, for example. Cohesive Zone Elements are typically defined at bonded contact interfaces between structural elements. The behavior of these Cohesive Zone Elements are governed by Traction Separation laws shown in Figure 1: as the cohesive zones separate, the traction first increases until a maximum interfacial allowable stress is reached, then the traction decreases to zero at a maximum interfacial allowable displacement, which results in complete local separation. Normal (Mode I) and shear (Mode II) behavior are both controlled by this triangular Traction Separation law. The area under the triangular curve is the Strain Energy Release Rate, and this energy dissipation rate can be experimentally determined for various joints.

Figure 1. Traction Separation Laws of Cohesive Zone Elements We will use NEi Nastran for the finite element solver, taking advantage of its threedimensional layered solid element capability. Like standard layered shell elements, the layup/stackup for the laminate is defined by the materials, ply thicknesses, and fiber orientation angles (if applicable) for each layer. A known limitation of layered shell elements is the usage of beam shear theory to derive balancing transverse shear stresses 2


on a ply by ply basis. Layered solid elements do not have this limitation – transverse shear and the through-thickness normal stresses are included as a direct result of the triaxial stress state in the fully 3D orthotropic material. The through-thickness stresses are taken into account to assess fiber failures, matrix cracking, and delamination between the plies. These stress components are particularly important when the following conditions are present: thick laminates and joints, out-of-plane loading, or the free edges of a component are the main regions of concern. With the same number of elements, there is no computational time penalty for using layered solid elements instead of layered shell elements. Also, the layered solid elements have three integration points per layer in the through-thickness direction to ensure accuracy of results. We will utilize NEi Nastran’s embedded MultiContinuum Technology (MCT) capability, provided by Firehole Composites’ Helius product. The traditional FEM approach for fiber reinforced composites is to smear the properties of the fiber and matrix together to arrive at a homogeneous representation of each layer. The smeared properties can be found from testing or a micro-mechanical approach. Failure in smeared-layer analyses is determined using a failure index that relates the stresses and strains in this smeared layer to experimentally-determined material strengths. The micromechanics-based MCT separates the matrix and fiber properties, handling each composite constituent differently and allowing failure to progress through a multi-step damage mechanism. For example, compression of a composite often results in the matrix being squeezed around with deformations due to the Poisson’s effect that drag the fibers with them. With the matrix and fibers in different stress states, the composite laminate can be much closer to failure than a smeared approach would predict. Effects like this can become significant especially when thermal loads are added to structural ones. Failure theories like Puck [5] and LARC02 [1] account for the different failure states of the fibers and matrix. They predict a two-step failure, with matrix (or possibly fiber) failing first, followed later by the other with a subsequent degradation of properties in the interim. MCT accounts for multiple-step failure, as well as out-of-plane stresses. Out-of-plane stresses can sometimes be a significant fraction of in-plane stresses, up to the point that the matrix can fail in the out-of-plane direction before the in-plane directions have failed completely. MCT also has the major advantage of being developed for weaves/fabrics and other composites (metal matrix, long fiber, unaligned fiber, etc.), while LARC02 and PUCK was developed strictly for unidirectional composite materials. Fertig provides the equations used to implement MCT [3]. Progressive Ply Failure Analysis allows a composite to fail one ply at a time. Once an iterative nonlinear model is setup (either nonlinear static or nonlinear transient dynamic), a Progressive Ply Failure Analysis can be run to examine the progression of nonlinear damage. The selected Failure Theory predicts when a particular ply will fail, and then the stiffness of the failed directions will be reduced with the corresponding reduction scale factor (10 percent is used for this investigation). Multiple papers have shown Progressive Ply Failure Analysis results to not be very sensitive to the exact value of the stiffness reduction scale factors – varying the reduction scale factors from 4-20 percent generally affect the composite FEA results minimally. Equilibrium is reestablished after the material properties are degraded, and then at the next time step/load increment, composite 3


stresses are again checked to see what plies and what direction the plies fail in. After the nonlinear analysis is concluded, the analyst can postprocess output such as percentage of plies failed and location of plies failed.

2. DOUBLE CANTILEVERED BEAM MODE I FRACTURE TEST 2.1 Lamina Properties The laminate for the Double Cantilevered Beam (DCB) tests are prepared from a carbon (T300)/toughened epoxy (977-2) prepreg supplied by Cytec Fiberite. Table 1. Material Properties for T300/977-2 Carbon Fiber

2.2 Finite Element Model Details A DCB FEM is constructed to explore the delamination propagation through a bonded joint with a starter crack – see Figure 2. The carbon fiber laminate has a stackup of [0]24, and has dimensions of 150mm long, 20 mm wide, and 4 mm thick. The starter crack is 55 mm long. Enforced displacement in the +Z/-Z direction is applied to the top edge of the top beam and the bottom edge of the bottom beam to enable a Mode I opening type fracture. Using 6 CHEXA (linear hexahedral solid elements) elements through the thickness, a total of 3600 CHEXA elements form the mesh. The elements are not connected in any fashion at the initial precrack locations, and cohesive zone elements are specified at the bonded contact interface. For the triangular Cohesive Zone Law, the Strain Energy Release Rate is related by: (1) 4


SDMAXT is the interfacial allowable stress, and UDMAXT is the interfacial allowable displacement. Multiple combinations of SDMAXT and UDMAXT were iterated on to fulfill the requirement of SERR = 0.3 N/mm [2]. For most combinations of SDMAXT and UDMAXT, FEA results were similar and not very sensitive. SDMAXT of 20 MPa and UDMAXT of 0.03 mm were found to provide good nonlinear convergence and correlate to experimental data well. Lowering SDMAXT and increasing UDMAXT usually helps with nonlinear convergence. Lowering SDMAXT excessively will result in significantly altered results. Progressive damage (PPFA) of the carbon fibers and epoxy matrices are ignored for this simulation.

Figure 2. DCB Finite Element Model with Cohesive Zone Elements 2.3 DCB Mode I Fracture Results Nonlinear static analysis (Nastran SOL 106) is run with large displacement effects that include follower forces and differential stiffness. The laminate opens up in a pretty stable manner – see Figure 3. The red regions indicate locations where the cohesive zones have completely failed and the adhesive is no longer functional.

Figure 3. Deformed Configuration of DCB (1:1 Actual Def)

5


Figure 4. Results Comparison – Simulation vs. Experiment Delamination initiates at a certain critical load (P Crit), and as the delaminations propagate through the adhesive bondline, the stiffness of the DCB decreases. Looking at Figure 4, the simulation results have good correlation with the experimental and analytical results. The critical delamination failure load from experiment is 58 N, and the failure load from Node 153 in the simulation is 53.6 N. The experimental results are taken from De Morais, et a. [2], and is a typical load-displacement curve from the testing. A critical load (P Crit) was also calculated to be 63.7 N through analytical relations – these calculations are not presented here. Nodes 151-155 are the 5 loaded nodes on the edge of the top beam (153 is in the center of the edge, 152 and 154 are next to 153, 151 and 155 are on the ends of the edge). Nodes 151 and 155 have the same load-displacement response, while nodes 152 and 154 have the same load-displacement response. Load-displacement response is plotted for all 5 nodes to show the variance of behavior through the width of the specimen. Experimental testing has shown that at the crack tip in DCB specimens, the Strain Energy Release Rate decreases close to the edges. The maximum load for the ends of the edge (Nodes 151 and 155) is 39.7 N.

6


3. MIXED MODE SKIN STIFFENER TENSION TEST 3.1 Lamina Properties Both skin and flange are made from IM6/3501-6 graphite/epoxy prepreg tape with a nominal thickness of 0.188 mm. Along with lamina material properties, constituent (fiber and matrix) material properties are also input into the finite element simulation for use with the MCT micromechanical failure model. Table 2. Material Properties for IM6/3501-6 Carbon Fiber [4]

3.2 Finite Element Model Details A bonded skin and flange assembly FEM is constructed to explore the complex damage mechanism of a complete delamination of a composite flange – see Figure 5. The skin layup, consisting of 14 plies, is [0/45/90/-45/45/-45/0]S, and the flange layup, consisting of 10 plies, is [45/90/-45/0/90]S. The skin width is 25.4 mm, and the skin length is 127.0 mm. Maximum flange length at the adhesive bond is 50 mm, and the minimum flange length is 42 mm. The angle at the flange tip is 27 degrees.

7


One side of the composite skin is constrained in X, Y, and Z translations. The other side is constrained transversely (Y and Z directions) and enforced to displace in the +X direction. Using 2 CHEXA (linear hexahedral solid elements) elements through the thickness for both the skin and the flange, a total of 828 layered solid elements form the mesh. Cohesive Zone Elements are specified for the bonded contact interface between the flange and the skin. The two tapered flange ends are meshed differently to prevent simultaneous delamination from occurring at both ends. For Mode I opening mode, SDMAXT is set to 11.3 MPa and UDMAXT is set to 0.0162 mm to satisfy SERR = 0.0931 N/mm [4]. For Mode II shearing mode, SDMAXS is set to 11.3 MPa and UDMAXS is set to 0.084 mm to satisfy SERR = 0.467 N/mm. For these types of CZM simulations, experience suggests we use a realistic value for interfacial allowable stresses and then adjust, if needed, to improve convergence. Micromechanics based damage degradation is utilized in this simulation.

Figure 5. Bonded Skin and Flange FEM 3.3 Skin Stiffener Tension Test Results

3.3.1 Displacement-Load Behavior Nonlinear static analysis (Nastran SOL 106) is run with large displacement effects that include follower forces and differential stiffness. NEi Nastran’s layered solid element and MultiContinuum Theory (MCT) capabilities add virtually zero computational penalty to a typical nonlinear finite element run – the biggest factor for run time is model size (number of degrees of freedom) and number of cohesive zone elements. The system is setup for large complex multiscale composite models, with simultaneous in-plane and out-of-plane damage progression, to run within hours. This particular analysis ran through 392 load increments, 9 load bisections, and 1625 equilibrium iterations in about half an hour (depending on hardware used). The displacement-load curve stays linear until matrix cracking occurs at the flange tip, which quickly propagates into debond failures along the adhesive bondline. A load drop occurs as the flange completely and unstably separates from the composite skin. In Figure 6, the x-axis is the displacements in the X direction extracted at the flange tip, and the yaxis is the total load applied in the X direction at the enforced motion boundary.

8


Figure 6. Displacement-Load Behavior of Skin Stiffener FEM Five tension tests were performed by Krueger [4]; the FE simulation results are tabulated with the experimental data in Figure 7. The correlation is good - the damage initiation load and flange debond load is slightly lower for the FE simulation, but for both simulation and test, damage initiation occurs at or above 90 % of flange debond load. It is interesting to note that for one test specimen, no visible damage could be detected prior to flange debond at 21,100 N. Figure 7. Results Comparison – Simulation vs. Experiment

9


3.3.2 Bond Failures and Matrix Failures From Figure 8, it is observed that the specimen is elongated by the tensile load, but also bends due to the local effects of the flange. The locations of matrix cracking and branching, after the complete separation of flange from skin, are shown in Figure 10. The delaminated bondline regions are shown in Figure 9.

Figure 8. Deformed Configuration at Damage Initiation (10 to 1 Actual Def)

Figure 9. Delamination Failure Regions

Figure 10. Matrix Failure Regions

10


3.3.3 Assessing Probability of Fiber Damage From using the micromechanics-based failure model MCT for damage degradation, we can extract matrix-specific and fiber-specific failure indices for the laminated structure at any load level. For example, at damage initiation load, matrix failure indices and fiber failure indices for all the plies are shown in Figures 11 and 12. Using this type of constituent-based information, the analyst can conclude that at damage initiation load, the fibers are not close to being damaged with a max Fiber FI of 0.22. Also, with the flange tips cracking, most of the flange exhibits matrix failure indices less than 0.1. The matrix cracking that occurs is in the through-thickness direction, and the middle of the composite flange does not experience much load.

Figure 11. Matrix Failure Indices at Damage Initiation

11


Figure 12. Fiber Failure Indices at Damage Initiation

4. CONCLUSIONS Composite bonded joints historically have been difficult to analyze for structural integrity and damage tolerance, and current aerospace programs have needed to depend heavily on expensive experimental testing to catch design flaws. An integrated Finite Element Analysis solution in the NASTRAN environment has been proposed for aerospace engineers to efficiently and accurately optimize complex composite structures. The combination of Cohesive Zone Modeling with 3D micromechanics based damage degradation offers invaluable design insight into the complex damage mechanism of composite bonded joints. For the test cases described in this paper, DCB Mode I Fracture Test and Skin Stiffener Tension Test, NASTRAN has proved to be fairly accurate in predicting critical delamination propagation failure loads. The failure load for DCB Mode I Fracture Test is predicted within 10 % of experiment, and the failure load for Skin Stiffener Tension Test is predicted within 20 % of experiment. Future research will focus on using this same coupling of in-plane and delamination simulation to model advanced dynamic load events – Barely Visible Impact Damage (BVID) and high speed ballistic impacts.

12


5. REFERENCES 1. Dávila, C.G., Jaunky, N., and Goswami, S. “Failure Criteria for FRP Laminates in Plane Stress.” 44th AIAA/ASME/ASCE/AHS Structures, Structural Dynamics, and Materials Conference, Norfolk, Virginia, 7-10 April 2003. 2. De Morais, A.B., De Moura, M.F., Marques, A.T., De Castro, P.T., “Mode-I Interlaminar Fracture of Carbon/Epoxy Cross-Ply Composites”. Composites Science and Technology 2002;62:679-686. 3. Fertig, Ray S. III, “Bridging the Gap Between Physics And Large-Scale Structural Analysis: A Novel Method For Fatigue Life Prediction Of Composites” 2009 SAMPE Fall Technical Conference, Wichita, Kansas, 19-21 October 2009. 4. Krueger, R., M. K. Cvitkovich, T. K. O’Brien, and P. J. Minguet, “Testing and Analysis of Composite Skin/Stringer Debonding under Multi-Axial Loading,” Journal of Composite Materials, vol. 34, no. 15, pp. 1263-1300, 2000. 5. Puck, A. and Schürmann, H. “Failure analysis of FRP laminates by means of physically based phenomenological models.” Composites Science and Technology,” 62:12-13 (2002) 1633-1662.

13


Cohesive Zone Modeling

Recommend Documents