July 2013 

NASA/TM–2013-218024 

Analysis Methods for Progressive Damage of 

Composite Structures 

Cheryl A. Rose, Carlos G. Dávila, and Frank A. Leone 

Langley Research Center, Hampton, Virginia 
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

Access the NASA STI program home page 

at http://www.sti.nasa.gov

E-mail your question to help@sti.nasa.gov

Fax your question to the NASA STI 

Information Desk at 443-757-5803 

Phone the NASA STI Information Desk at 

443-757-5802 

Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320

 
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

July 2013 

NASA/TM–2013-218024 

Analysis Methods for Progressive Damage of 

Composite Structures 

Cheryl A. Rose, Carlos G. Dávila, and Frank A. Leone 

Langley Research Center, Hampton, Virginia 
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not 

constitute an official endorsement, either expressed or implied, of such products or manufacturers by the 

National Aeronautics and Space Administration. 
Contents 

1 Introduction 2

2 Modeling Multiple Damage Mechanisms and their Interactions 1

2.1 Cohesive Models for Fracture in Composite Structures 4

2.1.1 Modeling Delamination Propagation in the Presence of Fiber Bridging 5

2.1.2 Improved Cohesive Models for Mixed-Mode Fracture 6

2.1.3 Cohesive Laws for Adhesive Modeling 7

2.1.4 Cohesive Elements for Improved Computational Efficiency 8

2.1.5 Extensions Required for Fatigue Loading 8

2.1.6 Summary, Cohesive Models for Fracture in Composite Structures 8

2.2 Continuum Damage Mechanics Models 9

2.2.1 Material Property Characterization 10

2.2.2 Limitations of Continuum Damage Mechanics Models 15

2.2.3 Refinements to CDM Models 18

2.2.4 Summary, Continuum Damage Mechanics Models 20

2.3 Discrete Damage Modeling Methods 21

2.3.1 Mesh Sensitivity/Refined Cohesive Elements for Proper Interaction 23

2.3.2 Crack Initiation, Saturation and Delamination Onset 24

2.3.3 Extensions for Compression after Impact and Fatigue Loading 25

2.3.4 Summary, Discrete Damage Models 25

3 Experimental Methods 26

3.1 Measuring Initial Imperfections with Coordinate Measurement Machine 27

3.2 Ultrasonic Inspection 28

3.3 Digital Image Correlation 29

3.4 X-Ray Computed Tomography 30

3.5 Summary, Experimental Methods 31

4 Applications 32

4.1 Damage Tolerance of Post Buckled Stiffened Panels 32

4.2 Open-Hole Tension 34

5 Summary 37
2 

1 Introduction 

 This document provides an overview of recent accomplishments and lessons learned in the 

development of general progressive damage analysis methods for predicting the residual strength 

and life of composite structures. The developments described were accomplished by the authors 

and university collaborators while supporting the Subsonic Rotary Wing (SRW) Project, the 

Aircraft Aging and Durability (AAD) Project and the Vehicle Systems Safety Technologies 

(VSST) project. These developments are described within their State-of-the-Art (SoA) context 

and the associated technology barriers. The emphasis of the authors is on developing these 

analysis tools for application at the structural level. Hence, modeling of damage progression is 

undertaken at the mesoscale, where the plies of a laminate are represented as a homogenous 

orthotropic continuum. 

 There are several primary motivations for developing analysis methods to model onset and 

growth of damage in composites. One is to eliminate some of the testing requirements in the 

building block approach used for design and certification of composite primary structures. By 

replacing the tests of subcomponents and components normally performed during the 

intermediate stages of the building block process, significant savings can be realized in terms of 

costs and the time required to progress from preliminary design to service. The second motivation 

for developing these analyses is to enable parametric studies to assess the relative severity of 

possible damage modes in new design concepts. Such a predictive capability offers the potential 

for the development of more efficient concepts by enabling the exploration of the large design 

space of laminated composites, which is not possible with the current largely empirical approach 

for structural design and substantiation. An additional motivation is to provide tools for assessing 

the safety of as-built structures and the structural integrity in the presence of manufacturing 

imperfections or flaws and in-service damage. 

 To predict reliably composite structural residual strength and life, general purpose 

progressive analysis tools capable of representing all of the composite damage modes (matrix 

cracking, delamination, fiber kinking, fiber breakage, etc.) and their interactions are needed. Of 

particular relevance for structural applications is the development of an improved understanding 

of the significance of interactions between matrix cracking and delaminations on structural life 

and strength, and the development of analysis methods to represent these interactions. For 

instance, correct representations of the effects of matrix cracks and delamination on fiber kinking 

failure phenomena are essential for accurate prediction of compression after impact strength. In 

addition, several studies have demonstrated that matrix shear nonlinearity can influence these 

interactions, and must be included in simulations. Due to the complexity of the failure processes 

that must be captured by the predictive models, the implementation of these capabilities into 

robust and accurate analysis software can only be achieved in incremental steps. 

The following sections provide a brief examination of the major advances in the SoA of 

progressive damage analysis for composite structures. The aim of the present effort is to establish 

the ranges of validity of available models, to identify technology barriers, and to establish the 

foundations of the future investigation efforts. Such are the necessary steps towards accurate and 

robust simulations that can replace some of the expensive and time-consuming “building block” 

tests that are currently required for the design and certification of aerospace structures. 
1 

2 Modeling Multiple Damage Mechanisms and their Interactions 

Modeling progressive failure in composite materials and structures remains a challenge after 

years of extensive research due to the large number of damage mechanisms that must be 

considered and their complex interactions. Fracture in a composite structure is the result of the 

evolution and interactions of discrete damage events such as fiber/matrix debonding, matrix 

cracking, delamination between plies, and fiber failure. Additionally, the effect of shear 

deformation on crack onset has been observed experimentally [1] and analyses have found that 

such nonlinearities must be represented by the models for an accurate prediction of the evolution 

of interacting matrix and delamination cracks [2-4]. These damage modes evolve in various 

combinations that depend on the laminate geometry and stacking sequence and cause 

redistributions of stresses in the failing composite. Some combinations of damage may reduce 

local stress concentrations, while others may precipitate a structural collapse. Therefore, an 

analysis methodology capable of predicting structural strength and fatigue life must take into 

account damage initiation and propagation. 

However, the details of the mechanisms that lead to failure are not fully understood due to the 

complexity involved with the idealization of the individual constituent responses and their 

interactions. The presence of two constituents, fiber and matrix, and the extreme anisotropy in 

both stiffness and strength properties result in damage mechanisms that act at different scales. 

The damage mechanisms can be divided into intralaminar and interlaminar damage. As shown in 

Fig. 1, intralaminar damage mechanisms correspond to fiber fracture and matrix cracking, 

whereas interlaminar damage mechanisms correspond to the interfacial separation of the plies 

(delamination). 

 a) Fiber fracture b) Matrix transverse crack c) Delamination 

Fig. 1. Damage mechanisms in laminated composites. 

The formulation of the governing physical principles of damage evolution depends on the scale 

of the idealization of the damage process, which may span from molecular dynamics to structural 

mechanics, and includes the intermediate scales of micro- and meso-mechanics. The damage 

models shown in Fig. 2 illustrate four typical scales of damage idealization. The micromechanical 

model shown in Fig. 2a represents what is normally the lowest practical scale of composite 

damage idealization, in which detailed matrix damage mechanisms, such as matrix plasticity and 

cracking, as well as fiber/matrix interface cracking are typically represented using representative 

volume element (RVE) [5], or unit cell models. The representation of damage at this level is 

typically based on a reduction of the material stiffness. Hence, fracture is represented as a band of 

localized volumetric stiffness reduction, referred-to as “soft discontinuity” as opposed to a “hard 

discontinuity” where voids are represented by displacement discontinuities in the model. Due to 

computational constraints, micromechanical models are typically two-dimensional and they 

represent domains much smaller than a ply thickness. Consequently, they are useful for 

representing a composite hardening response, before cracks localize at the ply level or larger 

scale. 
2 

Fig. 2. Levels of damage idealization, from micromechanical to structural. 

The most common idealizations of composite damage for representing structural level response 

are at the mesoscale level, where the individual plies are represented as homogenous, anisotropic 

materials (Fig. 2b and Fig. 2c). At this level, the intralaminar and interlaminar damages processes 

have traditionally been treated using different theories. Intra-ply damage modes have been 

investigated primarily within the framework of Continuum Damage Mechanics (CDM) [6], where 

damage modes are represented as a reduction in the corresponding stiffnesses. Delamination has 

been studied extensively using interface fracture modeling techniques such as cohesive zone 

models [7] and virtual crack closure techniques (VCCT) [8, 9]. Direct coupling of the intralaminar damage models with the interlaminar damage models is accomplished by incorporating 

both damage models within a single finite element model. Despite significant advances in 

progressive damage modeling, recent studies [10, 11] have indicated that CDM models for ply 

failure, coupled with cohesive zone models for delamination have difficulty in accurately 

representing laminate failure sequences that are characterized by strong coupling between matrix 

cracking and delamination. Among the limitations of the CDM models that contribute to this 

difficulty, are their inability to accurately predict matrix crack paths, and to describe accurately 

local effects of stress redistribution in a damaged area. 

To address the limitations in the CDM framework, recent emphasis has been placed on 

development of discrete damage models (DDM) at the meso-scale level with discrete 

representation of cracks of all types, so that their coupling can be directly accounted for in a 

prediction. Progress in this area has been enabled by nonlinear cohesive zone models of fracture, 

and critical developments in computational finite elements, based upon the eXtended Finite 

Structural

(Hard discontinuity) • Fracture mechanics and 

modifications

• Strain softening, cohesive 

laws

a.

Mesoscale

(Soft discontinuity) • Continuum damage 

models (CDM)

b. Mesoscale

(Hard discontinuity) • xFEM + Cohesive laws

c.

Micromechanical

(Soft discontinuity) • RVE models (unit cell)

d.

Scale of idealization Damage type Typical approaches

Through-the-thickness 

crack or delamination

Intralaminar damage

Discrete damage

Fiber/matrix interface, 

matrix plasticity and 

damage
3 

Element Method (X-FEM), that allow arbitrary cracks to initiate and propagate within a 

simulation, independently of the finite element mesh. 

Finally, the model of a through-crack in a fuselage panel shown in Fig. 2d illustrates an 

example of a structural-level damage model. The crack is represented as a hard discontinuity, and 

prediction of the propagation of the crack could be based on a strain-softening law or a criterion 

based on the critical energy release rate [12]. However, any structural-level crack propagation 

criterion is strongly dependent on the material system and laminate configuration and, 

consequently, must be determined for each new material system and laminate stacking sequence. 

In addition, structural-level semi-empirical fracture models cannot address the characteristics of 

the crack-tip damage zone nor the complex interactions between micro- and macro-failures 

associated with the crack-extension process. Instead, the crack-tip damage zone is simulated as 

some “effective” notch-tip damage zone that is assumed to grow in a self-similar manner. In 

many cases, self-similar crack growth is not observed, and the lack of resolution in the damage 

mechanisms often renders structural-level models inaccurate after a short propagation of damage. 

It is clear from the preceding overview of typical damage modeling strategies that the 

conceptual idealization of damage, i.e., the identification, characterization, and formulation of the 

governing physical mechanisms that constitute damage evolution, are different at each scale of 

idealization. Damage idealizations with higher resolution and kinematic freedom can capture 

multiple damage mechanisms with a separate damage law for each mechanism. These damage 

laws are likely to be simpler and require a fracture toughness that is lower than required in the 

damage laws used in coarser models to represent the same global response. For instance, mesoscale discrete-damage models may only need one simple cohesive law to represent a variety of 

matrix damage patterns, while CDM models need multiple empirical stiffness degradation laws 

and interacting activation functions to represent the softening of the material. 

Two computational schemes are being pursued for developing the analysis capability required 

for predicting strength and life of composite structures. In one approach, a CDM model for inplane damage is coupled with cohesive elements for interface damage. The second approach 

employs a DDM method for matrix cracking, in combination with cohesive elements for 

delamination. A critical distinction between the CDM methodology and DDM interface fracture 

models exists in the approach in which a displacement discontinuity is represented; i.e. the CDM 

methodology replaces the displacement discontinuity with local volumetric stiffness degradation, 

or a “soft discontinuity,” whereas the interface fracture based techniques directly include the 

kinematics of the displacement jump. The “soft discontinuity” approach used in the CDM models 

leads to some limitations of the CDM models, as described above. However, there are significant 

advantages of the CDM methodology which include relative maturity, general formulations and 

implementations in commercial finite element codes, the ability to describe in a computationally 

tractable manor the stiffness response of a laminate with large amounts of matrix damage, and the 

ability to represent all failure modes, including the onset and evolution of fiber failure which does 

not follow a predefined path. 

Recent applications of the DDM models [2, 4, 13, 14] have demonstrated high fidelity in their 

predictions for situations in which CDM approaches are limited, specifically failure scenarios 

dominated by matrix and delamination cracking. However, demonstrations to date have been for 

monotonic tension loading only; extension and validation for other loading cases requires 

extensive testing. Additionally, the developments for composite structures where multiple cracks 

in individual layers interact with each other, with multiple cracks in adjacent layers, and with 

delaminations between layers, are relatively immature. Furthermore, typical implementations are 

two-dimensional and do not include bending deformations, implementations in finite commercial 

finite element codes are limited, and extensions required to include interactions with fiber failure 

are not clear. Finally, the fine meshes required for representing the nonlinear fracture process 

zones of multiple cracks and their interactions may render these models intractable for 

simulations beyond the coupon level. 
4 

The remainder of the document is organized as follows. First, cohesive laws for representing 

nonlinear fracture processes are briefly described. These cohesive laws are used in the 

formulations of the cohesive elements for interface fracture, and for predicting the initiation and 

propagation of matrix cracks in the DDM models. In the DDM models, the same cohesive laws 

are used for matrix and delamination cracking. Recent advancements in determining complex 

cohesive laws for fracture processes where multiple energy dissipating mechanisms contribute to 

crack growth resistance are described, along with limitations for predicting crack growth under 

mixed-mode crack conditions. Extensions required for addressing fatigue loading, and for 

incorporating compression loading are also described. Then, the concepts of the continuum 

representation of composite material response are discussed and the intrinsic limitations of CDM 

models for laminated composites are outlined. A summary of recent developments to address 

these limitations and developments in experimental methods for determining critical material 

property input are presented as well. Next, an emerging DDM modeling technique based on the 

eXtended Finite Element Method (X-FEM) that overcomes many of the limitations of CDM 

models is presented. Capabilities of the CDM and DDM modeling techniques are then illustrated 

with some examples. Finally, advancements in experimental methods that have contributed to the 

validation of the methods are described, along with some suggestions for future required 

developments. To address the overall objective of the document, each section closes with a 

bulleted summary of the SoA, recent developments, and known technology barriers relevant to 

that section. 

2.1 Cohesive Models for Fracture in Composite Structures 

The ability to predict crack propagation in composites emerged four decades ago with the 

development of computational methods based on the theory of linear elastic fracture mechanics 

(LEFM). However, LEFM is limited to applications in which the fracture process zone is 

confined to the immediate neighborhood of the crack tip itself, or the fracture process zone is 

small compared to the characteristic crack length. For many fracture processes in composite 

materials and structures, the fracture process zone may be relatively large compared to other 

structural dimensions. The fracture process zone is a nonlinear zone characterized by plastic 

deformations and progressive material softening due to nonlinear material deformations, such as 

microcracking, void formation and fiber/matrix pullout. The size of the fracture process zone is 

dependent on the type of material softening and it must be considered in many situations of crack 

growth in composite structures. For example, development of a process zone gives rise to stable 

growth and crack growth resistance. The apparent fracture toughness increases with crack growth 

– an effect called the R-curve – until the process zone is fully developed. 

 Nonlinear fracture mechanics (NLFM) provides a framework for characterizing crack growth 

resistance and for analyzing initial amounts of stable crack growth. The cohesive zone model is a 

NLFM methodology that was developed to simulate the nonlinear fracture response near the 

crack tip, by collapsing the effect of the nonlinear process zone onto a surface of displacement 

discontinuity. Additionally, cohesive crack models have the ability to describe the process of 

void nucleation from inclusions. Therefore, they can be applied to initially un-cracked structures 

and can describe the entire fracture process, from no crack to complete structural failure. 

Cohesive crack models are based on kinematic descriptions that use discontinuities in the 

displacement field. Cohesive interfacial laws are phenomenological mechanical relations between 

the tractions and interfacial separations such that, with increasing interfacial separation, the 

tractions across the interface reach a maximum, and then decrease and vanish when complete 

decohesion occurs. It can be shown by performing a J-integral calculation along a contour 

surrounding a notch tip that the resulting work of interfacial separation is related to Griffith’s 

fracture criterion [15]: 
5 

 (1) 

where Gc is the critical energy release rate, and  and  are the interfacial traction and interfacial 

separation, respectively. 

Cohesive zone models have been used extensively in composite fracture analysis to represent 

multiple cohesive mechanisms acting across delamination or other types of cracks. For many 

common fracture processes in composite materials and structures, the apparent fracture toughness 

increases with crack growth. Such a response is typically due to the presence of more than one 

physical phenomenon involved in the separation process: some acting at small opening 

displacements, which are confined to correspondingly small distances from the crack tip, and 

others acting at higher displacements and extending further into the crack wake. In the presence 

of an R-curve, the toughness measured during crack propagation typically increases 

monotonically until reaching a steady-state value. In the case of delamination, the increase in 

toughness with crack growth is attributed to fiber bridging across the delamination plane. 

2.1.1 Modeling Delamination Propagation in the Presence of Fiber Bridging 

Fig. 3. R-curves for DCB specimen obtained by cohesive element superposition. 

The interlaminar crack growth resistance of long fiber-reinforced composites can significantly 

increase in the presence of fiber bridging, which induces tractions that act over an extended zone 

in the wake of the crack tip. As a consequence, the R-curve of unidirectional composites in Mode 

I delamination is characterized by a fracture toughness that gradually increases from an initial 

value of G1 up to a value corresponding to steady-state propagation, Gc [16, 17], as shown in Fig. 

3. To avoid the difficulties associated with R-curve effects, standard test methods such as ASTM 

D5528 [18] conservatively recommend the use of fracture initiation values to characterize the 

toughness of double cantilever beam (DCB) specimens. However, the use of the energy release 

rate measured at the initiation of nonlinearity causes a severe underestimation of the predicted 

strength, regardless of whether linear elastic fracture mechanics or cohesive zone models are 

used. Conversely, the use of the critical energy release rate for steady state propagation leads to a 

severe overestimation of the strength [16]. Therefore, for accurate predictions of strength, it is 

necessary to have a model that accounts for the R-curve effect. 

   



G d f

c   0
6 

The complex fracture processes involved in the delamination of composites can be 

approximated by superposing two or more different cohesive laws. The use of linear and nonlinear cohesive laws to model fiber bridging in delamination has been explored by many authors. 

It has been shown that softening laws within cohesive elements can be used to represent bridging 

as well as the dominant fracture process. The superposition of two bilinear cohesive responses 

was applied to model fiber bridging in the tensile fracture of unidirectional composites in Dávila 

et al. [16]. The trilinear cohesive law that results has also been used as an approximation of the 

complex fracture processes involved in the delamination of composites [19-22]. The principal 

difficulty in the use of superposed cohesive laws to represent R-curve effects is the lack of 

validated methodologies for the calibration of the softening laws. In addition, the development of 

a procedure for material parameter identification from DCB test results is complicated by several 

factors. It is known that the R-curves obtained are dependent on the specimen geometry and that 

there are no direct procedures for extracting bridging laws from experimental R-curves. Finally, 

the extraction of energy release rate (ERR) using standard procedures may be inaccurate for 

problems with large-scale bridging. 

In a recent study [23], experimental measurements were conducted to determine the critical 

ERR of glass/epoxy and its associated R-curve effect using DCB specimens. Then, the use of trilinear cohesive laws obtained by the superposition of simpler bilinear cohesive laws for 

representing the R-curve was investigated. Two approaches for determining the parameters of the 

superposed cohesive laws were proposed. In the first approach a procedure for extracting the 

cohesive parameters from experimental R-curves through the use of a new semi-analytical 

equation was developed. The methodology was applied to the experimental DCB data and it was 

shown that the resulting model was able to reproduce well the force-displacement response as 

well as the desired R-curve. In the second approach, an experimental R-curve is not required, and 

the parameters were obtained using a numerical optimization procedure that reduces the error 

between the predicted and experimental force-displacement results. The second approach is 

advantageous when fiber bridging introduces inaccuracies in the experimental energy release rate 

measurements. In addition, the second approach can be extended to identify the parameters 

needed for using more complex approximations of the cohesive laws. 

2.1.2 Improved Cohesive Models for Mixed-Mode Fracture 

In structural applications of composites, crack growth is likely to occur under mixed-mode 

loading. The loads and the number of cycles necessary to propagate constrained cracks such as 

delaminations are strongly dependent on the mode mixity, i.e. the ratio of the opening and 

shearing energies used in the formation of a new crack surface. Therefore a general formulation 

of cohesive laws must address mixed-mode fracture. Despite the maturity of cohesive laws, some 

issues regarding the prediction of crack propagation under mixed-mode conditions remain 

unresolved. It has recently been observed that although the global fracture mode mixity is 

constant throughout loading, the local fracture mode mixity is almost never constant during the 

damage evolution. Experiments by Högberg et al. [24] and analyses using cohesive elements 

show that the instantaneous mode mixity changes during the process of damage. It was found that 

damage initiation is usually dominated by transverse shear (Mode II) and, as damage evolves, the 

fracture process becomes increasingly dominated by Mode I opening. Consequently, the history 

of the interfacial traction, 	°, versus displacement, ∆ƒ

, relationship is not bilinear, as intended, but 

rather transitions from one bilinear function to the next, resulting in a concave curve as shown in 

red in Fig. 4. Unfortunately, all cohesive formulations assume that the mode mixity is constant 

during damage evolution. Therefore, cohesive and fracture models may not predict the correct 

propagation of fracture. 

In specimens such as the Mixed Mode Bending (MMB) specimen or the Double Cantilever 

Beam specimen with Uneven Bending Moments (DCB-UBM), where the fracture mode mixity 
7 

from a linear elastic viewpoint is independent of crack length, it is observed that all points along 

the crack path undergo the same history of mode change during damage. 

Fig. 4. Predicted instantaneous and global traction displacement laws in local mixed 

mode. 

An investigation is ongoing to understand better the significance of the dependence of the 

fracture propagation on the history of the mode mixity. The initial approach to this investigation 

consisted in a mathematical evaluation of the conditions under which cohesive element 

formulations are thermodynamically consistent. It was found that by imposing the so-called 

“Turon Constraint” [25, 26] a mixed-mode cohesive model can be derived from an energetic 

potential, therefore satisfying the laws of damage irreversibility and energy dissipation. However, 

such a model does not always tend to the linear elastic solution, even in cases of short process 

zones, as initially expected. It was also found that mixed-mode cohesive laws can reproduce 

exactly the analytical results only in specific cases which result in path independence. Additional 

work is needed to understand which results are correct and what cohesive models ensure that 

crack propagation is correctly predicted for all loading conditions and materials. 

2.1.3 Cohesive Laws for Adhesive Modeling 

A methodology using Digital Image Correlation (DIC) was developed to characterize the 

mixed-mode fracture properties of composite joints bonded with Cytec FM-300M mat-reinforced 

structural adhesive [27]. A test campaign was conducted using three types of specimens. Double 

Cantilever Beam (DCB), End Notch Flexure (ENF) and Mixed-Mode Bending (MMB) tests were 

performed on carbon-epoxy specimens to determine the Mode I, Mode II, and mixed-mode 

fracture properties, respectively. 

The characterization procedure uses DIC to extract the displacement histories around the initial 

crack tip, including the rigid body rotations as well as the shearing and opening displacement 

jumps across the crack tip. Following an inverse methodology, these displacement jumps are 

substituted into analytical J-integral equations that describe the fracture energy in terms of the 

rotations and displacement jumps at the crack tip. Finally, the cohesive laws needed for finite 

element analysis are obtained by numerical differentiation of the J-integral functions. The 

resulting cohesive laws represent the adhesive stiffness, its strength, and its fracture resistance in 

the form of stress-displacement functions in each mode of fracture. 

The prediction of mixed-mode fracture of adhesives is particularly challenging because it is 

difficult to satisfy the “Turon Constraint” [26] that ensures the correct energy dissipation. It was 
8 

found that without the imposition of this constraint, that the mode mixity is improperly computed 

by the elements and the predicted results are inaccurate. Additional research is necessary to 

develop a mixed-mode cohesive model that is not restricted by the assumptions of the Turon 

model. 

The measured cohesive laws were approximated using piecewise linear functions for numerical 

implementation into User defined cohesive ELements (UEL). The results from the numerical 

simulations of the test specimens confirm that the measured cohesive laws in combination with 

the UEL can be used to predict the debonding process in mixed mode, including the initiation of 

damage as well as steady-state tearing. 

2.1.4 Cohesive Elements for Improved Computational Efficiency 

One of the major drawbacks of cohesive elements is their requirement for fine meshes. Recent 

investigations have attempted to quantify a priori the mesh requirements and to establish 

procedures to determine the parameters that describe the cohesive laws. However, many 

structural problems remain computationally intractable due to these mesh requirements. Research 

is needed to identify solutions to reduce these mesh requirements. The solutions might require 

hybrid cohesive/VCCT formulations, and/or cohesive elements with advanced integration 

schemes. 

2.1.5 Extensions Required for Fatigue Loading 

Delamination can occur due to the accumulation of damage under both static and cyclic loads. 

Constitutive models for decohesion must be able to describe both types of material degradation. 

Models that extend the cohesive zone modeling approach for monotonic loading into forms 

suitable for cyclic loading, such as Turon [28], have not yet been thoroughly validated. In 

addition, very few models have been proposed that are truly predictive, i.e., that depend on 

measurable properties such as the Paris Law coefficients rather than relying on parameter 

calibrations that may be dependent on mesh refinement, mode mixity, R-ratio, Gmax, etc. Much 

work is needed to characterize experimentally the phenomena that lead to fatigue crack 

propagation, and to develop phenomenological laws that capture these effects. 

Finally, cohesive models always require fine meshes. However, in the case of fatigue loading, 

the mesh requirements quickly become intractable. Additional work is needed to reduce the 

requirements to have such excessively refined meshes. As in the case of quasi-static analyses, 

solutions that would alleviate the mesh requirements include hybrid cohesive/VCCT 

formulations, and cohesive elements with advanced integration schemes. 

2.1.6 Summary, Cohesive Models for Fracture in Composite Structures 

2.1 Cohesive Models for Fracture in Composite Structures 

State of the Art

Cohesive element technology has good degree of maturity as a 

robust prediction tool. Cohesive elements are becoming available 

in most commercial FEM software suites. 

Recent Advances

Procedures have been developed to determine a priori the 

parameters that describe the cohesive laws and the maximum size 

of the finite elements needed for proper energy dissipation. 
9 

2.2 Continuum Damage Mechanics Models 

The most common idealizations for representing the initiation and evolution of intralaminar 

damage modes in structural level predictions rely on meso-scale level continuum damage 

mechanics (CDM) models. At the meso-scale level, individual plies of a laminate are assumed to 

be composed of a homogeneous material with orthotropic properties. Continuum damage 

mechanics models for composite materials were pioneered by Ladevèze, Matzenmiller [31, 32] 

and others based on previous work by Kachanov, Lemaître [33, 34] and others. In-house analysis 

efforts have used primarily the Progressive Damage Model for Composites provided in Abaqus 

Relationships between the shape of a cohesive law and the R-curve 

effects are better understood. 

Progress made towards thermodynamically-consistent cohesive 

law formulation: definition of “Turon Constraint.”

Measurement of cohesive laws for adhesive failure in mode I and 

mode II demonstrated. 

Technology Barriers

Mixed-mode fracture poses challenges that are not well 

understood. Issues of path dependence must be addressed. 

Fatigue cohesive formulations are immature. 

Fatigue cohesive models have mesh requirements that are 

intractable. 

LaRC Cohesive Element Research Software 

Coh_UEL.03 (LaRC, [7]) 

· 8-node cohesive element with 3 DOF/node based on LaRC 

B-K criterion for mixed-mode fracture. 

· Secant stiffness only. 

· User-written UEL subroutine for Abaqus/Std. 

Coh_UEL.06 (LaRC, Turon, U. of Girona, [29]) 

· Thermodynamically-consistent constitutive model. 

· Analytically-derived tangent stiffness. 

Coh_UEL_Fatigue.06 (Turon, U. of Girona, [28]) 

· Fatigue damage model based on Turon’s PhD 

dissertation. 

Coh_UEL_Shell.07 (LaRC, [30]). 

· 8-node cohesive element with 6 DOF/node allows 

connection of shell elements without requiring coincident 

nodes. 

Coh_UEL.12 (Sarrado/Turon U. of Girona, [26]) 

· Numerical tangent stiffness for ease of constitutive model 

development and improved convergence rate. 

· Interlaminar shear stiffness is calculated internally for 

thermodynamic consistency in mixed-mode. 
10 

[35] and the continuum damage model originally proposed by Maimí [36]. In these composite 

damage models, a distinction is made between the different failure modes, especially between 

fiber and matrix failure, and damage is represented by progressively degrading appropriate 

stiffnesses in accordance with predefined failure criteria. There are many failure criteria that have 

been proposed to predict the onset of matrix cracking and fiber fracture, some of which are 

described and compared in the World Wide Failure Exercise [37, 38]. However, few criteria can 

represent several relevant aspects of the failure process of laminated composites, e.g. the increase 

of apparent shear strength when applying moderate values of transverse compression, or the 

detrimental effect of the in-plane shear stresses in failure by fiber kinking. The LaRC03 failure 

criteria [39] and subsequent evolutions [40] address some of the limitations of other failure 

criteria as identified from the WWFE. For example, the LaRC criteria account for the effect of 

ply thickness, fiber misalignment in compression, and the effect of shear nonlinearity on fiber 

kinking and in-situ strength. 

In the model proposed by Maimí [36], the LaRC04 failure criteria [40] are used as damage 

activation functions, FM, to formulate a continuum damage model to predict the propagation of 

the several (M) damage mechanisms occurring at the intralaminar level. Each damage activation 

function predicts one type of damage mechanism using the following equations: 

 (2) 

where are internal variables (equal to 1 at time t = 0 ), and the functions correspond to 

the LaRC04 failure criteria. When a damage activation function is satisfied, FM =0, the associated 

damage variable, dM, is greater than zero, and the ply compliance tensor is affected by the 

presence of damage. Once a damage initiation criterion is satisfied further loading causes 

reduction of the material stiffness coefficients. The reduction of the stiffness coefficients is 

controlled by the damage variables, whose evolutions in the post-damage initiation phase are 

represented by softening laws.

In many damage models, including the Progressive Damage Model for Composites provided in 

Abaqus [35] and most cohesive laws, the evolution of damage is represented with bilinear laws 

that are described by a maximum traction and a critical energy release rate. The computational 

implementation of softening constitutive equations poses difficulties, however, because the 

boundary value problem becomes ill-posed when strain localizes along a surface known as the 

failure surface, leading to mesh subjective schemes [41]. Upon reducing the element sizes to zero, 

these analysis models predict failure to occur with zero energy dissipation. To resolve this lack of 

objectivity with respect to element size, a characteristic length is introduced into the constitutive 

model using a procedure based on the crack band model proposed by Bažant and Oh [42]. A 

complete definition of a continuum damage model for the simulation of intralaminar damage can 

be found in [43]. 

2.2.1 Material Property Characterization 

The Abaqus damage model and the damage model described in the 2006 NASA TM [36] use 

as input ply-level material properties and fracture toughness values. The use of ply-level 

properties is advantageous because standard test methods are available for determining the 

majority of the material properties required for input to the damage model, imperfections and 

variabilities that must be considered at lower scales can be disregarded, and use of ply-level 

properties eliminates testing requirements to determine material property input for the damage 

models every time the lay-up or stacking sequence is changed. These models require as input the 

material properties below: 

:   
 0 t

M

t

M M F   r

t rM  t M 
11 

Ply elastic properties (E1, E2, G12, G23, v12, v23) and ply strengths (XT, XC, YC, SL). 

Four components of the fracture toughness associated with longitudinal failure in tension 

and compression (G1+ and G1-, respectively) and transverse cracking in mode I and II, 

(G2+ and G6, respectively) 

All of the ply elastic and strength properties, with the exception of the ply shear strength, can 

be determined using test standards defined by the American Society for Testing Materials 

(ASTM). Traditionally, due to limited characterization data, it has been assumed that predictions 

of the growth of intralaminar cracks can be based on the fracture toughness values for growth of 

interlaminar cracks. The mode I fracture toughness component for transverse matrix cracking, 

G2+, is therefore determined using the Double Cantilever Beam (DCB) test ASTM-D 5528, and 

G6 can be determined using the four-point bending end-notched flexure test [44]. Recent 

experimental studies to determine mode I intralaminar fracture toughness from compact tension 

tests have confirmed that the intralaminar and interlaminar mode I fracture toughness values are 

equal [45]. 

 Resistance Curves in Fiber Failure of Composites. There are no standard test methods for 

determining the fracture toughness associated with longitudinal failure in tension and 

compression. It has been recommended that the fracture toughness for tensile longitudinal failure 

be determined using the compact tension test developed by Pinho [46]. Recently a new 

methodology was proposed to measure the crack growth resistance curves (R-curves) associated 

with fiber-dominated failure modes directly from these tests. The proposed methodology is based 

on identification of the crack tip location using the Digital Image Correlation (DIC) technique and 

the calculation of the J-integral directly from the test data using a simple closed form expression 

derived explicitly for the specimen lamination sequence [47]. The method obviates the need of 

any complex pre- and post-processing of the test data, either based on finite element calculations 

or standard data reduction methods, and enables the real-time generation of R-curves during a 

test. 

The crack-growth resistance curves obtained not only characterize the fracture toughness of the 

material, but they also provide the basis for identification of the parameters of the softening laws 

used in the numerical simulation of the longitudinal fracture in the continuum damage mechanics 

models. The shape of the softening law is often assumed to be inconsequential for prediction of 

fracture, provided that it is defined as a function of the fracture toughness. While this assumption 

is valid when crack propagation is governed primarily by a single energy dissipating mechanism, 

the shape of the cohesive law plays a critical role in the prediction of fracture when crack 

propagation includes multiple energy dissipation mechanisms that act over different length scales. 

In this case, each energy-dissipating mechanism must be appropriately accounted for in the 

cohesive law. 

Several failure mechanisms, including fiber tensile fracture, fiber-matrix pullout and matrix 

cracking are present when a crack propagates in a plane perpendicular to the fiber direction. It 

was demonstrated by Dávila et al. [16] that a simple bilinear law is unable to predict the load 

displacement response obtained in a cross-ply compact tension test. A tri-linear law was 

proposed, with parameters of the softening law determined from the experimentally determined 

R-Curve, using the procedure defined by Dávila, et al. The load-displacement response observed 

in the test was successfully predicted when the determined softening law was used in a simulation 

that represented the specimen fracture with cohesive elements embedded in the model along a 

fracture plane extending from the initial notch tip.

Predictions obtained with the softening law incorporated in the continuum damage mechanics 

model, however, were not as successful. Preliminary simulation results, shown in Fig. 5 predict 

the development of splitting cracks parallel to the 0-degree fibers fairly early in the loading. The 
12 

cracks develop over large portions of the coupon, and are not constrained to the immediate 

vicinity of the idealized fracture surface. These cracks likely developed in the test as well, 

however, their presence is not explicitly noted, and their contribution to the measured R-curves is 

not clear. The possible presence of these splitting cracks in the compact tension specimens affects 

both the total energy dissipated into the creation of new fracture surfaces and the severity of the 

stress concentration ahead of the notch. Furthermore, the measured fiber fracture toughness 

becomes dependent on the matrix properties and stacking sequence of the cross-ply specimens. 

Fig. 5. Matrix splitting predicted via CDM in a compact tension finite element model.

Results of these simulations illustrate the difficulty in properly interpreting the experimental 

data, and backing out the necessary input for the damage models. As more and more detail and 

fidelity in representing evolving damage mechanisms are incorporated into the damage models, 

care must be taken to prevent double-counting the contributions of individual energy dissipating 

mechanisms. Experimental studies with detailed NDE examination, including in-situ ultrasound, 

and X-Ray/CT during interrupted tests, to examine the evolution of all of the damage modes 

during the coupon loading, coupled with high-fidelity analysis are required to resolve these 

issues. To this end, several compact tension specimens, with varying stacking sequences would 

be required to decouple the measured fiber fracture toughness from the specimen configuration 

dependencies. New NDE techniques with resolution sufficient to capture fiber failures and fibermatrix debonding may be required as well to provide sufficient characterization of all of the 

significant damage modes and their interactions. 

Determining the fracture toughness for longitudinal failure in compression is more challenging. 

Compact compression tests were previously proposed to be a good configuration for measuring 

resistance curves associated with fiber kinking. However, delamination that accompanies the 

propagation of a kink-band renders the compact compression test unsuitable for this purpose. Test 

specimen and experimental design are required to determine the fiber fracture toughness in 

compression. 

In-Situ Strengths and Shear Nonlinearity. To predict matrix cracking in a laminate subjected 

to in-plane shear and transverse tensile stresses, a failure criterion must account for the ‘in-situ’ 

strengths. Within multidirectional laminated composites, laminar strength oftentimes does not 

correspond directly to the experimentally measured strengths of a unidirectional composite. 

Adjoining plies of differing orientations constrain the laminae, affecting their effective strength. 

These in-situ strengths vary with the unidirectional strength, fracture toughness, and thickness of 
13 

the individual plies. The in-situ effect, originally detected in Parvizi’s [48] tensile tests of crossply glass fiber reinforced plastics, is characterized by higher transverse tensile and shear strengths 

of a ply in a laminate when it is constrained by plies with different fiber orientations, compared 

with the strength of the same ply in a unidirectional laminate. The in-situ strength also depends 

on the number of plies clustered together, and on the fiber orientations of the constraining plies. 

The results of Wang’s [49] tests of [0/90n/0] carbon–epoxy laminates indicate that thinner plies 

exhibit higher transverse tensile strength.

The in-situ strengths have been solved for by Camanho and Dávila for both matrix tensile and 

shear modes [50]. Due to the different fracture mechanisms acting on plies of different 

thicknesses and at different positions within a laminate, different in-situ strength equations were 

developed for thin embedded, thick embedded and surface plies. For the case of thin embedded 

plies, the general equation used to solve for the in-situ strength is: 

ܩୡ ൌ గ௧

౪ౢఌ౫ ߝߪ ׬ ସ

଴ (3) 

where ܩୡ is the fracture toughness, ݐ is the ply thickness, and ߝ୳୪୲ is the strain corresponding to 

the unidirectional strength of the material. Using equation (3), all strain energy preceding failure 

(both elastic and inelastic) is related to the corresponding experimental fracture toughness of the 

material through a fracture mechanics analysis of a constrained ply. In-situ strengths are found to 

have an inverse relationship with ply thickness (or the number of plies blocked together), 

increasing with decreasing thickness. The in-situ strength for thick embedded plies, however, is 

constant, not scaling with thickness. 

Several continuum damage mechanics models have been proposed that either neglect shear 

nonlinearity or provide the end-user with the ability to ignore the effects of shear nonlinearity in 

their particular model. The simplification of neglecting shear nonlinearity, however, is at odds 

with the use of in-situ shear strengths. Shear nonlinearity is integral to the current derivation of 

in-situ shear strengths [50], and electing to neglect shear nonlinearity while utilizing in-situ shear 

strengths is contradictory. 

For example, as currently derived, in-situ shear strengths are highly dependent on the shape of 

the shear stress-strain response (e.g., linear vs. nonlinear) prior to failure. To avoid adopting 

contradictory assumptions, it is necessary to model the nonlinear shear behavior in the manner 

that is assumed during the derivation of the in-situ shear strengths. Furthermore, much effort has 

been devoted to properly evolving the state of matrix damage throughout the full range of 

possible mode mixities (see below). Even moderate amounts of shear nonlinearity can have a 

considerable influence on the calculated mode mixity for a given strain history, causing drastic 

changes in the effective toughness of the local material and the overall strength of the modeled 

structure. 

In their original derivation of in-situ shear strengths, Camanho and Dávila utilized a shear 

nonlinearity equation by Hahn and Tsai in which a single parameter is used to fit experimental 

shear data [51]. However, this equation was found by Schuecker, et al. [52] to lack the necessary 

flexibility to closely match all experimental shear responses. In order to maintain flexibility in the 

shape of the shear response curve, the two-parameter Ramberg-Osgood law [53] was selected to 

represent the nonlinear shear response: 

ߛଵଶ ൌ ଵ

ீభమ

ൣ߬ଵଶ ൅ ߙ୮୪ሺ߬ଵଶሻȁ߬ଵଶȁ

௡൧ (4) 

where ߙ୮୪ and ݊ are used to fit experimental shear data. A comparison of the abilities of the 

Hahn-Tsai and Ramberg-Osgood laws to fit experimental shear data from the second World-Wide 

Failure Exercise is shown in Fig. 6. Clearly, the added flexibility of the Ramberg-Osgood law 

allows for a much better approximation of the experimental data. 
14 

Fig. 6. Comparison of the Hahn-Tsai and Ramberg-Osgood shear nonlinearity 

equations with experimental shear data from the second WWFE. For the HahnTsai curve, ࢼ is equal to 2.80E−8; for the Ramberg-Osgood curve, ܔܘࢻ and ࢔

are equal to 4.41E−10 and 5.93, respectively.

The use of the Ramberg-Osgood shear law in place of the Hahn-Tsai law required that the insitu shear strengths be re-derived. The primary difference between the Ramberg-Osgood and 

Hahn-Tsai laws is the replacement of the cubic exponential term with the variable exponent ݊. 

The inclusion of this a priori unknown exponent ݊ causes there to be no closed-form solutions for 

the in-situ shear strengths. The Newton-Raphson method is used to obtain an approximate 

solution for the in-situ shear strength ܵ୐ǡ୧ୱ by solving equations (5) through (7) when the 

equations are equal to zero: 

Thin Inner Plies: ݂ሺܵ୐ǡ୧ୱሻൌ

Ͷܩଵଶܩ୍୍ୡ

ݐߨ െ ܵ୐ǡ୧ୱ

ଶ

ʹ െ ߙ୮୪݊

݊൅ͳ ܵ୐ǡ୧ୱ

௡ାଵ (5) 

Thick Inner Plies: ݂ሺܵ୐ǡ୧ୱሻ ൌ ቆ

ܵ୐ǡ୧ୱ

ଶ

ʹ ൅

݊୮୪ߙ

݊൅ͳ ܵ୐ǡ୧ୱ

௡ାଵቇ െ ʹ ቆ

ܵ୐

ଶ

ʹ ൅

݊୮୪ߙ

݊൅ͳܵ୐

௡ାଵቇ (6) 

Surface Plies: ݂ሺܵ୐ǡ୧ୱሻൌ

ୡ୍୍ܩଵଶܩʹ

ݐߨ െ ܵ୐ǡ୧ୱ

ଶ

ʹ െ ߙ୮୪݊

݊൅ͳ ܵ୐ǡ୧ୱ

௡ାଵ (7) 

An example of the predicted variation of the in-situ shear strengths with ply thickness for an 

embedded ply of IM7/8552 is shown in Fig. 7, using three different pre-failure shear responses: 

linear, Hahn-Tsai, and Ramberg-Osgood. Clearly, the in-situ shear strengths are highly sensitive 

to the posited shape of the shear response, with the in-situ strength decreasing with increasingly 

nonlinear responses. Also worth noting is the much lower transition between the thick and thin 

ply regions: for the case of the Ramberg-Osgood law fit to the WWFE data, a single 0.125-mm 

thick ply of IM7/8552 is considered thick, and, therefore, the strength is not predicted to scale 

with ply blocking. 

0

20

40

60

80

100

0.0 1.0 2.0 3.0 4.0 5.0

Shear Stress [MPa]

Shear strain [%]

Test (WWFE)

Hahn-Tsai

Ramberg-Osgood

Material: IM7/8552
15 

Fig. 7. In-situ shear strengths for IM7/8552, as a function of ply thickness. The 

corresponding unidirectional shear strength is 92.3 MPa. 

For a material with a given fracture toughness, the low in-situ shear strengths corresponding to 

highly nonlinear material responses lead to very high in-situ ultimate strains. For example, a 

single embedded 0.125-mm thick ply of IM7/8552 will have in-situ ultimate strains of 5.5%, 

8.9% and 9.8% using the linear-elastic, Hahn-Tsai, and Ramberg-Osgood shear laws, 

respectively. For the case of thin-ply laminated composites (e.g., 0.04-mm thickness) [54], these 

strains increase to 9.7%, 19.7%, and 24.6%. These strain values are far beyond the experimental 

unidirectional results (approximately 5% for IM7/8552). It is expected that extrapolations of the 

nonlinear shear response curves are invalid at these very high strains. Without further 

experimental data demonstrating the relationships between nonlinear shear behavior and in-situ 

shear strengths in various laminate configurations, the reliability of predicted in-situ shear 

strengths is questionable. 

The significant dependence of in-situ shear strengths on the assumptions made regarding the 

shape of the nonlinear shear stress-strain response highlights the importance of maintaining 

consistency between the selection of in-situ shear strengths and the assumed shear behavior for a 

particular model. Furthermore, a greater understanding of the nonlinear shear response of 

IM7/8552 is required, including, but not limited to, its dependence on the in-situ effect and on 

transverse loading. 

2.2.2 Limitations of Continuum Damage Mechanics Models 

Despite advances in progressive damage modeling, recent studies [10, 11] indicate that CDM 

models coupled with cohesive zone models may not always represent laminate failure sequences 

properly. These deficiencies are particularly evident when the observed fracture mode exhibits 

matrix splitting and pullouts [55] or when the fracture is characterized by a strong coupling 

between transverse matrix cracking and delamination. The deficiencies of the predictive 

capabilities consist of several issues, including the incorrect prediction of the damage zone size 

normal to the fracture direction when using crack-band models and the inability of local CDM 

models to reliably predict matrix crack paths. These limitations are mostly due to the fact that 

CDM models are usually implemented as “local” rather than “non-local” models [56], i.e., the 

evolution of damage in a local CDM is evaluated at individual integration points without 

0

100

200

300

400

0 0.2 0.4 0.6 0.8 1

In-situ Shear Strength [MPa]

Ply Thickness [mm]

Linear

Hahn-Tsai

Ramberg-Osgood

UD
16 

consideration of the state of damage at neighboring locations. The following discussion pertains 

mostly to such local implementations, since non-local damage models are less widely used due to 

the difficulty in implementing them within the finite element method. 

The premise of the crack-band approach for regularizing CDM models is that damage localizes 

into a band with a width equivalent to the element dimension. If the element size is smaller than 

the damage process zone, the crack-band approach may not predict correctly the width of the 

damage zone nor the local stress field. Consequently, the stress redistribution resulting from 

damage development may be inaccurately predicted and can potentially result in inaccurate 

representation of damage mode interactions and failure sequences. 

As a result of homogenization and damage localization, CDM models have difficulty 

predicting crack paths. Since homogenization eliminates the distinction between fibers and 

matrix, a CDM model cannot distinguish between cracks that propagate along fiber directions 

from those that cross fibers when loaded in shear [10]. In CDM models implemented with 

damage localization, the damage state at any integration point in the model depends only on the 

stress field at that point rather than the damage state of neighboring points. Therefore, the 

direction of damage evolution is driven only by the instantaneous local stress distribution. In 

other words, the local direction of cracking may be predicted correctly by the failure criteria, but 

the sequence of failures that eventually defines the path of a crack at a macroscopic level may be 

predicted incorrectly. 

The inability of CDM models to determine the correct direction of propagation is particularly 

evident when the stress field is dominated by shear. Consider two different plies in a laminate 

with a notch and subjected to shear, as shown in Fig. 8 In both situations, the direction of the 

matrix microcracks is correctly predicted by the failure criterion to be 45°. Furthermore, both 

situations would result in an identical sequence of failures, since the stress field is identical. 

However, it is clearly not the same to propagate a crack in a sequence of linked microcracks (Fig. 

8a) as it is to propagate a crack band across fibers (Fig. 8b). Matrix cracking in a shear band 

running parallel to the fibers is a relatively brittle failure mechanism, whereas matrix cracking 

normal to the fibers produces a damage band that requires much more work to propagate. 

 a) Propagation of shear damage along fiber b) Propagation across fiber direction 

Fig. 8. Idealized propagation of shear damage (adapted from [10]). 

The sensitivity of CDM predictions to the finite element mesh orientation also contributes to 

the difficulty in predicting the crack path. Although the objectivity of the solution with respect to 

element size is addressed with the crack-band approach described previously, the predicted 

damage may be dependent on mesh orientation and element shape. When strain-softening 

constitutive models are used in a finite element simulation, damage tends to propagate along 

preferred directions, consisting of either element edges or element diagonals. A demonstration of 

the sensitivity of simulation results to mesh orientation is provided in Fig. 9 for a unidirectional 

compact tension specimen with the fiber placed at 90° to the load direction. Results are presented 

	12

Damage band

	12

a

	12

	12
17 

for simulations obtained with a mesh oriented parallel to the fiber direction (Fig. 9a) and with an 

inclined mesh in front of the crack tip (Fig. 9b). The crack should propagate along the fiber 

direction. However, the results show directional bias, and the simulated crack band propagates in 

the direction of the element alignment. 

 

 a) Mesh aligned with crack direction b) Mesh inclined to crack direction 

Fig. 9. Effect of mesh orientation on crack path in a unidirectional CT specimen. 

The tendency for damage to localize along mesh lines can be partially attributed to shear 

locking [56]. In the CDM methodology, a crack or displacement discontinuity is represented by a 

degradation of the corresponding terms in the constitutive stiffness. As the crack opens, the 

stiffness degradation is such that stress should not be transferred across the crack faces. However, 

such unloading may not occur due to in-plane shear locking. Shear locking here refers to 

inappropriate shear stress transfer across a widely open smeared crack, which occurs when an 

element cannot shear without inducing tensile strains. In a ply with orthotropic properties, a 

matrix crack should be represented by setting the transverse shear modulus, G12, and the 

transverse Young’s modulus, E22, to zero. However, quadrilateral elements have been shown to 

exhibit coupling between 
12 and 11 unless the element edges are aligned with the softening band 

or are oriented at 45 degrees to the band [56]. Furthermore, the tendency of the element to lock is 

dependent on the order of integration of the element: fully integrated elements are more 

susceptible to pathological in-plane shear locking than reduced-integration elements. The shear 

stress transfer across an open smeared cracked can result in inaccurate prediction of stress 

redistribution after damage development. 

Load transfer across CDM cracks is also related to the complexity of representing the 

kinematics of a fully or partially damaged material with a single integration point. CDM models 

are implemented at the integration point level. Typically, a single integration point is used to 

represent the constitutive response of a single material. After the local initiation of damage, 

however, an integration point is tasked with representing the evolution of a crack and the 

constitutive responses of the bulk material on each side of the crack. While the local stress terms 

that are normal to and acting on the new crack plane will be equal in each of these three material 

sections, the stress terms acting parallel to the crack plane in the remaining bulk material will be 

artificially linked due to their shared integration point. Ideally, the stress terms parallel to the 

CDM crack could be solved for and stored separately for the material on either side of the crack, 

similar to two solid elements with a cohesive element in the middle. 

When damage localizes and a fracture path is known a priori, a relatively simple approach to 

circumvent some of the limitations of CDM models noted above consists of aligning the mesh 

with the direction of fracture [57] to force a matrix crack to localize along the fiber direction. The 

benefit of mesh alignment was demonstrated by Song et al. [58] for a 0-degree unidirectional 
18 

open-hole tension specimen, loaded parallel to the fibers, and in quasi-isotropic open-hole tension 

specimens whose failures were dominated by development of matrix and delamination cracks. 

Predictions obtained with aligned meshes were significantly improved over those obtained with 

a traditional radial mesh, and demonstrated an ability to predict properly crack paths, and failure 

sequences, where predictions obtained with traditional radial meshes were unsuccessful. 

However, the lack of kinematic freedom within a simulated crack represented by a band of 

softened elements prevented total stress relaxation across the open crack, leading to some 

inaccuracy in calculations of stress redistributions. Additionally, use of structured meshes 

corresponding to the fiber orientation in each ply to describe properly matrix cracks and thus 

reduce mesh dependency, leads to non-coincident meshes through the laminate thickness. Tie 

constraints are therefore required to connect the individual plies in a finite element model. The 

use of tie constraints results in large run times and computational resource requirements, which 

may render these types of models intractable for large applications. 

2.2.3 Refinements to CDM Models 

In addition to developing modeling approaches for addressing some of the deficiencies in 

CDM approaches for representing matrix and delamination cracking interaction, significant effort 

has been devoted to refining the continuum damage mechanics (CDM) finite element (FE) 

subroutine (LaRC CompDam CDM software [59]) corresponding to the damage model originally 

proposed by Maimí [36]. 

 Refinements to Continuum Damage Mechanics Model. The first modification addresses 

deficiencies in the CDM model presented in [36], for predicting crack propagation. The CDM 

model used separate matrix tension (or compression) and shear damage state variables to affect 

the damaged material response, each of which were functions of only the elastic limit. That is, as 

the evaluated mixed-mode failure criterion increased in value, the evolution of the damage 

variables was independent of the local mode mixity. A new mixed-mode matrix damage 

evolution law has been developed, based on the cohesive zone damage model of Turon [25]. In 

the updated CDM model a single matrix tension (or compression) and shear damage variable is 

used, defined using mixed-mode properties corresponding to the ratio of the transverse and shear 

strains [59].

In addition to the new matrix damage evolution law, changes have been made to the stresses 

which drive damage evolution and to the shape of the nonlinear shear responses prior to the onset 

of damage that can be modeled. Within the CDM model, a combination of the LaRC03 [39] and 

LaRC04 [40] failure criteria are used to drive the development of the damage state variables. 

These failure criteria are evaluated in terms of their effective stress. These effective stresses are 

representative of the stress state in the intact cross-sectional area of the element, and can 

accurately drive damage from onset to complete failure while a single damage mode is present. 

As was noted by Schuecker [52], however, the subsequent initiation and development of 

additional damage variables can lead to erroneous results with the current definition of effective 

stresses. To address this issue, the development of each damage variable is defined in terms of all 

other damage variables. For example, the effective stresses that drive the formation of the fiber 

damage variables consider the current matrix damage state, and the effective stresses that drive 

the formation of the matrix damage variables consider the current fiber damage state. This 

revision of the effective stresses removes the potential for spurious additional damage modes to 

be predicted and allows for a more accurate simulation of multiple simultaneous local damage 

mechanisms. 

The presence of any plasticity and/or damage at an integration point can lead to very large 

strains in a still partially intact volume of material. In order to properly represent the deformed 

material and the orientation of any local fracture surfaces while under finite strains, the LaRC 
19 

CompDam CDM software has been updated to perform all stress and strain calculations in terms 

of the 2nd Piola-Kirchhoff stress tensor and the Green-Lagrange strain tensor. 

Smeared Crack Formulation for Representing Matrix Cracks. An alternate formulation of 

the CDM model has also been developed that includes an embedded cohesive law for the 

representation of matrix damage within a continuum finite element, based on the smeared crack 

formulation of Camanho, et al. [60]. The goal of this implementation is to represent accurately 

the formation, opening, and closure of matrix cracks in fiber-reinforced composite materials 

under large deformations by improving the kinematics of cracked material. If matrix cracks are 

represented solely by softening certain terms of the constitutive stiffness matrix of damaged 

elements, the determined orientation of cracks within damaged elements becomes increasingly 

erred when the elements undergo large shear deformations. An incorrectly determined crack 

normal causes elements containing matrix damage to “collapse” normal to the fracture plane 

when under shear load, and contributes to spurious load transfer across the fracture plane due to 

the fiber direction being incorrectly influenced by the elemental shear. The inclusion of the 

embedded cohesive law provides improved bulk material stress and strain predictions, as well as 

improved predictions of crack opening or closure displacements under large shear deformations. 

To accomplish this goal, it is necessary to : i) establish the orientation of the crack; ii) develop 

a consistent means of describing the bulk and cohesive deformations in terms of the elemental 

deformation; and iii) determine the distribution of elemental deformation between the bulk and 

cohesive materials. The implementation is carried out by decomposing the elemental deformation 

gradient tensor into two parts: a new deformation gradient tensor associated with only the bulk 

material, and a displacement vector describing the cohesive displacement-jump. The orientation 

of the cracks is defined according to the bulk material deformation gradient tensor. The 

distribution of deformation between the bulk and cohesive domains is determined by considering 

the current stiffnesses of the two materials. Stress equivalence on the rotated fracture plane is 

assumed, and the Newton-Raphson method is used to solve for the distribution of bulk and 

cohesive displacements whenever there is a predicted change in the damage state. 

An example of the improved predictive capability offered by the smeared crack formulation is 

shown in Fig. 10, for a unidirectional open-hole tension specimen. While the accurate prediction 

of matrix cracks (Fig. 10a) is possible with traditional CDM methods, problems arise after 

continued stretching of the specimen, including nonphysical contraction normal to the fracture 

plane and increasing load transfer across the fracture plane with continued shear deformation. 

Using the cohesive-embedded CDM model, there is no increasing component of stress transfer 

across the matrix splits (Fig. 10b and Fig. 10c), and the deformation normal to the matrix fracture 

plane is maintained (Fig. 10d). 
20 

Fig. 10. Open-hole tension finite element model results: (a) matrix damage; (b) fiber 

stress; (c) constant stress concentration factor with increasing elemental shear 

deformation (~1.08); and (d) cracked element deformation under high shear 

deformation. 

2.2.4 Summary, Continuum Damage Mechanics Models 

2.2 Continuum Damage Mechanics Models 

State of the Art

CDM models that use failure criteria and regularized softening 

laws can predict the initiation and propagation of intralaminar 

damage modes. 

CDM is an improvement over Progressive Failure Analysis, which 

is mesh-dependent. 

Recent Advances

Progress made in identifying and addressing the limitations of 

CDM. 

Improvements made to LaRC CompDam CDM software, include: 

use of 2nd Piola-Kirchhoff stress tensor and deformation gradient; 

embedded cohesive crack within bulk continuum, new definition of 

effective stresses; improvements in mixed-mode matrix cracking 

prediction. 
21 

2.3 Discrete Damage Modeling Methods 

The eXtended Finite Element Method (X-FEM) is a technique that can be used to predict the 

location and evolution of matrix cracks in composites while avoiding the aforementioned 

limitations associated with continuum damage mechanics models. X-FEM is a mesh enrichment 

technique based on a pioneering concept by Moës [61] that facilitates the introduction of 

displacement discontinuities such as cracks at locations and along directions that are independent 

of the underlying finite element mesh. Although most of the research on X-FEM is devoted to 

arbitrary crack propagation in isotropic materials, recent applications to composite materials 

include delamination modeling and textile composite architecture representation [62, 63]. Huynh 

[64] provides a review of contemporary development of X-FEM, as well as novel applications to 

interfacial cracking analysis in two- and three-dimensional problems. 

Technology Barriers

Experimental characterization must be capable of separating 

different toughening mechanisms such as notch blunting by fiber 

bridging from fiber fracture according to the kinematic complexity 

of the model. 

Convergence rate of nonlinear solution in implicit analyses 

LaRC CDM Research Software

CompDam_UMAT.06 (LaRC, P. Maimí, U. Girona, [6, 43]) 

· Continuum Damage Model for laminated composites 

based on the LaRC03 damage activation functions [39].

· Two-part linear/exponential softening for fiber damage, 

exponential softening for all other damage modes. 

· User-written UMAT subroutine for Abaqus/Std. 

CompDam_UMAT.10 (LaRC, P. Maimí, U. Girona, unpublished) 

· Bilinear softening laws of fiber tension and fiber 

compression. 

· Linear softening laws for matrix tension and matrix 

compression. 

· Simple plasticity model in inplane shear 

· New multi-mode damage variable for postprocessing. 

CompDam_VUMAT.13 (LaRC, unpublished) 

· Damage-mode-dependent effective stresses for proper 

evolution of multiple simultaneous failure processes. 

· Finite strain and stress formulation used for representing 

wide-open cracks. 

· Mixed-mode matrix crack softening. 

· Improved shear plasticity modeling. 

· Transverse shear plasticity. 

· Implemented as a VUMAT subroutine for Abaqus/Explicit. 
22 

Fig. 11. Matrix crack and delamination interactions in a composite laminate: (a) initial 

stage without damage, (b) matrix cracking stage, (c) delamination stage, linking 

up of matrix cracks in various plies, (d) specimen fracture [72]. 

Modeling a matrix crack in a ply that propagates parallel to the ply fiber direction is 

conceptually straightforward using X-FEM. However, it is more difficult to model networks of 

matrix cracks in a laminate where the fracture planes of matrix cracks in individual plies intersect 

at common interfaces and can cause delaminations that link the cracks through the thickness, as 

shown in Fig. 11. Within the traditional X-FEM approach, the difficulty in modeling networks of 

linked matrix cracks could be addressed by developing a special enrichment for multiple crack 

situations or by connecting two enriched/cracked elements. Such connection was recently 

accomplished in a quasi-two dimensional formulation by van der Meer et al. [65] and Ling et al. 

[66]. 

Another direction in which the X-FEM is being developed is the regularized extended finite 

element method (Rx-FEM) proposed by Iarve [67-70]. In the Rx-FEM approach, the usual step 

function used in X-FEM approaches to describe the crack surface is replaced by a continuous 

function. Displacement shape functions are used to approximate the step function, and the Gauss 

integration can be retained for element stiffness matrix computation regardless of the orientation 

of the crack. The cohesive connection between two plies in which matrix cracks have been 

introduced can then be easily established by computing integrals of the products of the shape 

functions at the ply interface. Therefore, a kinematically powerful model of crack networks can 

be constructed in which transverse matrix cracks parallel to the fiber direction are inserted with 

the Rx-FEM technique at locations determined using a failure criterion. A cohesive interface 

damage model is used to represent matrix crack propagation and delamination between plies [71].

Preliminary validation studies of X-FEM methodologies for un-notched and notched laminates 

subjected to tensile loading have demonstrated high fidelity in predicting experimentally 

observed responses, indicating the potential of this methodology [2-4, 71, 72]. An example of the 

predictive capability of the Rx-FEM approach developed by Iarve et al., is shown in Fig. 12 for a 

laminated composite Overheight Compact Tension (OCT) specimen [2]. Predicted damage 

patterns shown on the left are compared with measured X-Ray/CT data [73] on the right for a 

[452/902/-452/02]s specimen loaded to a pin-opening displacement (POD) of 2.55 mm. The 

damage patterns shown in the figures are a superposition of all cracks (black in Fig. 12a, and 

white in Fig. 12b) and all delaminations (blue, green and orange) through the specimen thickness. 
23 

A comparison of the predictions and experimental results shows very good correlation between 

the simulated and experimental delamination shapes and extent, and the matrix cracking 

distribution.

2.3.1 Mesh Sensitivity/Refined Cohesive Elements for Proper Interaction 

Adopting a modeling approach as described above with the ply material modeled with 

advanced X-FEM elements to allow for arbitrary cracking, and cohesive elements at ply 

interfaces to account for potential delaminations, should ideally provide a direct coupling 

between developing ply cracks and interface cracks. However, when a discontinuity is inserted in 

the displacement field of a plane element that is connected to a cohesive interface element, the 

relative displacement between the planes is affected. Allowing a solid element representing a ply 

to crack, theoretically requires the interface element connected to the cracked ply to be modified 

accordingly to ensure proper stress transfer across the interface. Additionally, the possibility that 

cracking has occurred in elements above and below the interface plane should be accounted for. 

To address this issue, an augmented cohesive zone (ACZ) element has recently been proposed 

that is based on the augmented finite element (A-FEM [66]) formulation to allow for arbitrary 

separation of the cohesive element in accordance with the separations that have occurred in 

adjacent solid elements [74]. The effectiveness of the proposed ACZ elements has been 

demonstrated for several sample problems. However, the additional complexity in the numerical 

models may not be required. Standard non-augmented interface cohesive elements may be 

sufficient for most applications, especially when using fine meshes [71, 75]. Additional studies 

are required to better understand the modeling requirements and to validate modeling 

assumptions. 

 a) Predicted (b) X/Ray-CT 

Fig. 12. Superposed predicted and experimental cracking and delamination extent for a 

[452/902/-45/02]s compact tension specimen loaded to a POD of approximately 

2.55 mm. (From Mollenhauer [2], with permission) 
24 

2.3.2 Crack Initiation, Saturation and Delamination Onset 

Failure analyses with X-FEM are most often performed at the meso-scale level, in which each 

ply is represented by a homogeneous orthotropic material. To capture the proper sequence of load 

redistributions that result from interactions between matrix cracks and delaminations, each ply 

must be modeled with separate elements. Therefore, the minimum mesh requirement consists of a 

single layer of elements across the thickness of each ply and a layer of cohesive elements between 

the plies. Few analyses have been attempted with more than one element over the thickness of the 

ply because the computational requirements associated with such small elements render 

intractable any analysis larger than a small coupon. However, the elliptical opening profile of a 

transverse crack cannot be represented with a single element, as illustrated in Fig. 13, for a crossply [0/90n/0] laminate, loaded in tension parallel to the 0-degree fibers. Studies have been 

initiated to investigate the ability of meso-scale level X-FEM models with a single layer of 

elements per ply to capture accurately all aspects of matrix cracking. In particular, whether the 

models can predict the in-situ ply thickness effect on crack initiation and propagation, the crack 

density as a function of strain, the strain for crack saturation, and the interaction between 

delamination and transverse cracks are being assessed. 

In collaboration with F. van der Meer, TU Delft, it has been shown for the first time that 

initiation and propagation of transverse matrix cracks in composite laminates can be captured 

with cohesive zone analysis [76]. Experimental results show that the load levels at which the 

cracks grow in the transverse direction increase with decreasing thickness of the ply, a trend 

which is correctly predicted by the cohesive model. It was shown that a second critical load level 

exists which is associated with propagation in the thickness direction. For thin plies, the in-situ 

strength is governed by transverse propagation, while through-thickness propagation is critical for 

thick plies. 

Stress relaxation due to crack opening plays an important role in the transverse crack 

propagation, which governs the thin ply in-situ strength. This crack opening cannot be captured 

when only one element is used through the thickness of a ply and the interface between plies is 

rigid prior to delamination onset. However, a technique based on shear lag was introduced which 

makes the interface deformable. Crucial for a proper prediction of the in-situ effect is that the 

stiffness of the interface is inversely proportional to the thickness of the transverse ply. The 

required interfacial properties were determined analytically using a closed-form shear lag model. 

As a result, the dependence of the propagation load level on the thickness can be captured with a 

cohesive method with a single element per ply. Preliminary results indicate that this modified 

model can also capture the density of transverse matrix cracks as a function of the applied strain. 

Additional work is necessary to ensure that the interfacial compliances do not adversely affect the 

softening response of the interfacial cohesive elements. 

 (a) Multiple elements (b) Rigid interface (c) Deformable Interface 

Fig. 13. Side view of final deformation for three different cases: multiple elements 

across the thickness (a), single element with rigid interface (b) and single 

element with deformable interface (c); shading indicates displacement in load 

direction. 
25 

2.3.3 Extensions for Compression after Impact and Fatigue Loading 

The initial development of the X-FEM methodologies focused on modeling transverse matrix 

crack and delamination interaction for tension loading. Fiber failures were not considered. 

Preliminary validation studies for un-notched and notched laminates subjected to tensile loading 

have demonstrated high fidelity in predicting experimentally observed responses, indicating the 

potential of this methodology [2-4, 71]. NASA currently has a grant with the University of 

Dayton Research Institute (UDRI) to extend the methodology to include tensile fiber damage, 

matrix cracks that are inclined through the thickness of a ply, as develop in certain combinations 

of shear and compression loading, and fatigue loading conditions. Additional extensions required 

to address the problem of predicting damage due to impact include accounting for geometrically 

non-linear deformations, and transient loading conditions. Extensions to address the prediction of 

compression after impact strength include, as a minimum, coupling with a CDM model to include 

the effects of fiber on the compression failure. 

Additional development of the fatigue constitutive model described above is required to ensure 

that the same model can predict transverse matrix cracks within the X-FEM framework as well as 

delamination propagation in cohesive element models. For fatigue loading, material 

characterization testing will be required and may include additional Paris law coefficients, S-N 

curves, etc. Additionally, pathfinder experiments are required for both compression loading and 

fatigue loading to guide the model development, to determine the critical phenomena that must be 

modeled and at what scale. 

2.3.4 Summary, Discrete Damage Models 

2.3 Discrete Damage Models 

State of the Art

X-FEM modeling techniques are gaining wider acceptance. 

Commercial FEM vendors such as DS-Simulia/Abaqus and ANSYS 

are developing X-FEM capabilities. 

Current commercial capability usually allows a single crack in a 

domain, which is insufficient for composite analysis. 

Majority of formulations are two-dimensional; bending 

deformations are not included. 

Predictive capability demonstrated for limited applications 

subjected to monotonic tension loading. 

Recent Advances

LaRC has spurred the development of two X-FEM modeling 

techniques (A-FEM and Rx-FEM) that are capable of predicting 

complex crack networks in laminated composites and the failure 

that results from the interactions of the cracks and delaminations. 

Technology Barriers

Methods are needed to predict matrix cracking and the associated 

in-situ effects with a single element per ply. 
26 

3 Experimental Methods 

Recent advances in damage mechanics have led to increasingly complex models of damage 

evolution in composite structures with vast numbers of new computational models capable of 

predicting the damage processes in composites continuously being advocated. However, major 

challenges remain to improve the quantitatively predictive power available to the engineer today 

through these unprecedented computational/analytical techniques. The process of improving 

predictive capability requires a close interaction between experimentalist and theories to ensure 

that the basic building blocks of the theories are well established. 

In the past five years there has been significant interaction with the Nondestructive Evaluation 

Sciences Branch (NESB) at NASA Langley Research Center (LaRC) to use existing NDE 

techniques and to develop new experimental capability within the structural and materials 

laboratories to allow for detailed pre- and post-test inspection of test specimens, as well as 

inspection under load. Routine examinations prior to all tests include accurate measurement of 

surface geometry using coordinate measurement machines, and ultrasonic inspection to provide 

an indication of specimen quality. These data are input as initial conditions into computational 

models. For small coupons X-Ray/CT examination is also being performed. In-situ techniques, 

including acoustic emission, digital image correlation techniques, and non-immersion ultrasonic 

inspection are also being employed to monitor the onset and propagation of damage, and to 

provide an indication for test interruption and more detailed examination by X-Ray computed 

tomography (CT). The following examples illustrate the use of new experimental methods to 

improve the understanding of the failure processes in composites, to improve the predictive 

capability of the models, and to validate the models. 

X-FEM models are not capable of predicting fiber failure. Hybrid 

X-FEM/CDM models are needed. 

Extensions are required for compression damage, and for fatigue 

loading. 

Computational requirements limit application to small specimens. 

Methods for coupling with continuum representations are 

required. 

LaRC DDM Research Software

BSAM (U. of Dayton Research Institute (UDRI), LaRC, [71]) 

· Regularized extended finite element model for laminated 

composited based on X-FEM 

· Mixed-mode bilinear cohesive laws for matrix crack and 

delamination propagation [29], LaRC03 failure criteria 

for matrix crack insertion [39], and CompDam CDM 

model [6, 43]. 

· Stand-alone analysis code and VTMS post-processor 

A-FEM_UEL.10 (U. of Miami, Teledyne Scientific, LaRC [66]) 

· Augmented-FEM approach based on X-FEM and extended 

to laminated composites. 

· User-written UEL for Abaqus/Std. 
27 

3.1 Measuring Initial Imperfections with Coordinate Measurement Machine 

The response and failure of thin structures subjected to compression loads is strongly affected 

by small initial imperfections and residual thermal deformations. These initial deformations can 

be measured using a coordinate measurement machine. This device uses a laser probe to generate 

high-precision spatial coordinates of millions of points on the surfaces of the specimens. The 

precision of each point is 10-6 mm, which allows an accurate analysis of microscopic 

deformations. In the photograph shown in Fig. 14, the CMM is used to characterize the 

deformation of a single stringer compression specimen before it is tested to failure. The measured 

shape of the skin (Fig. 15a), can be compared to the predicted residual thermal deformation (Fig. 

15b). These initial shapes are then incorporated into the analysis model to be used for predicting 

the nonlinear response of the structure. 

Fig. 14. Measurement of geometric imperfections using the Coordinate Measurement 

Machine with a laser probe. 

 

(a) Measured imperfection (b) Predicted residual deformation 

Fig. 15. Measured imperfection and predicted residual deformation after thermal cool 

down. 
28 

3.2 Ultrasonic Inspection 

New non-immersion high-definition ultrasonic measurement (UT) systems are being routinely 

applied to inspect specimens before tests. The UT system, which is shown below uses a focused 

immersion probe mounted in a captive water column with a rugged Nitrile membrane tip at the 

focal point. Performing inspections before the tests ensures that any defects are properly 

accounted for in data interpretation and in analysis simulations. For instance, it was found during 

a routine scan of the specimen shown below that an anomaly was present along the stiffener 

flange (Fig. 16). The specimen was nevertheless tested to failure and its collapse load was 

significantly lower than expected. Without the initial scan, the manufacturing defect would have 

gone unnoticed and the results of the test could have been misunderstood. 

 

Fig. 16. Pre-test ultrasonic inspection of experimental test articles. 

A non-immersion system has also been developed for in-situ application [77] so that damage 

initiation and progression can be monitored while a specimen is under load (Fig. 17). This 

capability has been used during the testing of specimens loaded monotonically and for specimens 

loaded in fatigue until failure. The system uses normal-incidence pulse-echo ultrasound, which is 

optimal for detecting delaminations, but is not as good for characterizing matrix cracking. 

Methods employing non-normal incidence ultrasound, or polar backscattering, are being 

developed to provide a rapid capability for a more complete characterization of damage, 

including measurement of matrix cracking [78]. 

Defect
29 

Fig. 17. Non-immersion in-situ ultrasound (UT) system. 

3.3 Digital Image Correlation 

Digital image correlation (DIC) is being used during testing to track the onset and growth of 

surface damage, and to develop methodologies for determining crack growth as a function of 

load. Methods for determining sub-surface damage from changes in surface strain profiles are 

also being investigated. 

In a recent application, a methodology using the commercial measuring system Vic-3D [79] 

was developed to characterize the fracture of composite joints bonded with Cytec FM 300-M 

adhesive [80]. A test campaign was conducted using three types of specimens. Double Cantilever 

Beam, End Notch Flexure and Mixed-Mode Bending tests were performed on bonded carbonepoxy specimens to determine the mode I, mode II, and mixed-mode fracture properties, 

respectively. The characterization procedure uses DIC to extract the displacement histories 

around the initial crack tip (Fig. 18a-c). Following an inverse methodology, these displacements 

are substituted into analytical J-integral equations that describe the fracture energy (ERR), Fig. 

18d. Finally, the cohesive laws needed for finite element analysis are obtained by numerical 

differentiation of the J-integral functions. The resulting cohesive laws represent the adhesive 

stiffness, strength, and fracture resistance in the form of stress-displacement functions in each 

mode of fracture (Fig. 18). 
30 

Fig. 18. Procedure to estimate the Mode II cohesive law: (a) displacement field around 

the crack-tip without external load; (b) displacement field around the crack-tip 

at maximum load; (c) output from the interrogation of the vertical line through 

the crack tip; (d) J-integral vs. displacement-jump; and (e) the experimental 

cohesive law. 

3.4 X-Ray Computed Tomography 

Recent advances in X-Ray/CT capability have provided the opportunity for detailed 

examination of damage development and evolution in the interior of composite structures. 

Resolution on the order of 10-20 μm is relatively common with commercially available sources, 

and can be obtained for specimens of laboratory scale with reasonable data acquisition and 

reconstruction times. The ability to obtain three-dimensional images of damage patterns provides 

the data necessary to tremendously increase understanding of the mechanisms of damage 

development and interaction, and to quantify predictive capability in terms of damage 

development as well as ultimate load. Additional advances in ultra-high resolution CT (~ 0.5 μm) 

are yielding images of damage mechanisms with spatial resolutions below the ply thickness [1, 

81]. Data at this scale may be required for improving understanding of the mechanisms involved 

in the initiation and early stages of damage evolution, the interactions of matrix cracks with fiber 

failures, and may be critical for refinement of fatigue damage models. 

X-Ray/CT examination has become a critical component in the validation process of the 

damage models. Nondestructively imaging internal damage, during interrupted tests, provides the 

opportunity to monitor progression of damage with load for comparison with the damage 

predictions. A sample image shown in Fig. 19 obtained from an open-hole tension specimen 

shows the ability to resolve both delaminations and matrix cracks. Data acquisition and 

reconstruction times, however limit application of X-Ray/CT examination to small specimens or 
31 

very small regions of specimens for ultra-high resolution inspection. New X-Ray sources and 

significantly reduced data acquisition and reconstruction times are required for in-situ highresolution inspection and inspection of large parts. 

Fig. 19. X-Ray/CT to monitor internal damage development. 

3.5 Summary, Experimental Methods 

3 Experimental Methods 

State of the Art

In-situ digital image correlation, ultrasonic inspection (UT) 

routinely performed for laboratory scale specimens. 

X-Ray/CT capability with resolution on the order of 5-20 μm, 

common with commercially available sources, can resolve matrix 

cracks and delaminations. Data acquisition and reconstruction 

times are reasonable for laboratory scale specimens. 

Recent Advances

The Nondestructive Evaluation Sciences Branch (NESB) and The 

Structural Mechanics and Concepts Branch (SMCB) at LaRC have 

spurred the development of in-situ UT inspection capability, and 

rapid UT inspection capability for large structural subcomponents. 

Emerging stereo microscope digital image correlation systems 

provide order of magnitude increase in resolution over 

conventional lenses, and provide potential ability to detect 

precursors to crack nucleation. 

Matrix crack

Delamination

Crack

Delamination
32 

4 Applications 

4.1 Damage Tolerance of Post Buckled Stiffened Panels 

Rotorcraft fuselages are usually thin, so they are susceptible to buckling at low strain levels. 

Major gains in structural efficiency can be achieved by exploiting the ability of stiffened 

composite structures to operate deep into the postbuckling regime. These potential weight savings 

are not yet fully achievable, mainly because of the difficulty in predicting the strength of 

composite structures, and the difficulty in determining the durability and damage tolerance of 

postbuckled structures. Consequently, there is increasing interest in the study of damage tolerant 

structures. Using the specially-developed Single Stringer Compression Specimen (SSCS) shown 

in Fig. 20, experimental and analytical studies have been conducted to investigate the effects of 

bond defects, geometric imperfections, and postbuckling response on the collapse of stiffened 

composite panels. 

Predicting the failure of stiffened panels in postbuckling is difficult for several reasons. One of 

the difficulties is associated with numerical convergence when the nonlinearities induced by 

propagating material failures are interacting with geometric non-linearity as a result of the large 

post-buckling deformations. An additional difficulty is the sensitivity of the panel response to 

initial geometric imperfections, because the postbuckling response paths and the associated 

failure modes can differ among nominally identical panels due to the presence of the initial 

geometric imperfections. A study of the effect of geometric imperfections on the predicted 

response was conducted by applying small perturbations in the skin normal direction to the 

locations of the nodes in the finite element model. The maximum perturbation was kept at 1/10 of 

the skin thickness. The distribution of the perturbations was calculated by linear combinations of 

the first three buckling modes. It was found from a large number of analyses that different 

imperfections resulted in one of the three different responses shown in Fig. 21. The differences in 

the postbuckling responses are particularly important in the presence of defects or pre-existing 

damage. For instance, a bond defect between skin and stringer at the location shown in Fig. 21 is 

more likely to propagate if the deformation is of Type A, where the skin pulls away from the 

stringer, as opposed to Type C, where the deflections of the skin are small, or Type B, where the 

skin deflects into the stringer. 

The collapse of single stringer specimens is sudden, making it difficult to determine the 

sequence of events that lead to the catastrophic failure. High-speed video recorded during the 

tests confirms the numerical predictions: the postbuckling deformations cause a local pull-off of 

the stringer flange from the skin, as shown in Fig. 22. The ensuing delamination grows under the 

stringer, which becomes unstable and it fails in a crippling mode. The digital image correlation 

Emerging ultra-high resolution CT (~0.5 μm) provides spatial 

resolution below the ply thickness. 

Technology Barriers

Data acquisition and reconstruction times limit application of XRay/CT examination to small specimens, or very small regions of 

specimens for ultra-high resolution inspection. New X-ray sources 

and significantly reduced data acquisition times are required for 

high resolution in-situ inspection and inspection of large parts. 

In-situ ultrasonic inspection is in general limited to 

characterization of delamination. Developments are required in 

array-based polar backscattering approaches for characterizing 

transverse matrix cracks and delaminations. 
33 

system records the critical deformations of the skin, which are shown in Fig. 23. A change in the 

buckling mode of the skin inside the stringer initiates a pull-off (Fig. 23a and Fig. 23b) and 

results in a detachment of the skin along the entire specimen (Fig. 23c). 

Fig. 20. The Single Stringer Compression Specimen was developed to investigate the 

damage tolerance of postbuckled stiffened structures. 

Fig. 21. Three postbuckling response types predicted using different initial geometric 

imperfections. 

Damage Tolerance of Postbuckled Composite Airframe

Airframes consist of thin skins that are stiffened 

longitudinally by stringers and circumferentially by 

frames.

The design of thin structures subjected to 

compression is challenging, especially for composites.

Simplified shell model with intralaminar damage and 

delamination:

Composite fuselage

Out-of-plane response

A B C D E

Linear 

deformation

First buckling

Initiation of matrix damage

Initiation of 

delamination

Collapse

Delamination

Crippling Cr

Hat

Skin

Cohesive elements
34 

Fig. 22. Failure sequence in SSCS specimen exhibiting Type A postbuckling response. 

 (a) (b) (c) 

Fig. 23. Type B surface deformation sequence measured using digital image correlation. 

4.2 Open-Hole Tension 

The open-hole tension specimen is one configuration being used to assess the performance of 

the analysis methods that are being developed. This specimen, although relatively simple in 

geometry and loading, represents a significant challenge to progressive damage analysis tools, 

because of the complexity of the damage evolution and interactions that occur prior to ultimate 

failure. The open-hole tension (OHT) specimen was recently studied, both experimentally and 

computationally [55, 82]. Specimens were fabricated from IM7/8552 graphite-epoxy and had a 

quasi-isotropic stacking sequence. Three-dimensional scaling, with ply-level and sub-laminate-
35 

level approaches to through-the-thickness scaling, was examined to determine size effects in 

laminate failure. 

An in-house experimental study is being conducted to examine further the evolution of damage 

and failure in these specimens, and to assess the modeling approaches that are being developed. 

Current tests are focused on monotonic tension loading. Specimens have been fabricated for 

subsequent testing under cyclic loading. Ply-level scaled specimens only are considered in the 

present study because of the damage modes that develop in these specimens. Ply-level scaled 

specimens were fabricated by blocking multiple plies with the same orientation together, thereby 

increasing the effective ply thickness. Two thicknesses are considered with the stacking sequence 

[45m/90m/-45m/0m]s and with m = [2, 4]. The plate width to hole diameter ratio, w/d is kept 

constant at a value of 5, and the length to hole diameter ratio, l/d, is equal to 20. 

The experimental data indicate that the failure progression and the mode of laminate failure are 

dependent on the laminate layup (stacking sequence and ply thicknesses) and that the scaling of 

strength with respect to size is dependent on the type of failure. Sub-critical damage develops at 

relatively low loads prior to ultimate failure and typically begins with matrix cracks in the offaxis (non-load bearing) plies and progresses through interconnecting delaminations in a complex 

three-dimensional pattern, as shown in Fig. 24. It is the accumulation of these damage modes that 

ultimately leads to failure of the specimen. The extent of the sub-critical damage development is 

dependent on the laminate geometry and significantly influences the laminate failure mode. 

Three failure modes have been identified, and are classified as brittle, pull-out and delamination 

failure modes [55]. The brittle failure mode is relatively easy to predict, the remaining are 

difficult. 

The two computational schemes described above, CDM for in-plane damage coupled with 

interface elements for interlaminar damage, and DDM for in-plane damage coupled with interface 

elements for interlaminar damage have been employed for predicting the response of the OHT 

specimens. Complete details of the modeling approaches and analysis predictions can be found in 

Swindeman [72], and Song et al. [58]. The analysis conducted by Swindeman used the Rx-FEM 

methodology developed by Iarve [71]. The analysis conducted by Song used the standard Abaqus 

CDM model, and a finite element model with the meshes in each ply aligned with the fiber 

direction to address the CDM limitations described above. Analyses were also conducted by Song 

using a radial mesh in all plies to demonstrate the inability to predict the complex damage 

evolution using traditional modeling approaches. 

 Overall load versus end-displacement predictions agreed well with experimentally obtained 

data, and the effects of scaling on failure load were well predicted for the delamination failure 

mode. In addition, predicted damage patterns qualitatively agree well with experimentally 

obtained images of damage development. An example comparison of predictions and 

experimental data for a specimen with a 1/8-inch diameter hole, just prior to specimen failure is 

provided in Fig. 24. In all figures, delamination damage at the interface indicated is shown, along 

with transverse matrix cracks in adjacent plies. In the CDM images, gray areas are delaminations, 

and red lines are cracks. In the DDM predictions, shaded areas are delaminations, and lines are 

predicted cracks. In the experimental images gray areas are delaminations and black lines are 

cracks. As shown, both computational methods are able to predict well the complex damage 

patterns experimentally observed. The CDM predictions show fewer matrix cracks than the DDM 

predictions, and in the experimentally observed damage. These additional matrix cracks are in 

some cases significant and can allow for link-up of damage through the thickness of the laminate, 

earlier than the case with fewer matrix cracks. Additionally matrix cracks in the different plies are 

loaded in different mixed-mode conditions, and the eventual interaction of the cracks with 

delaminations and cracks in adjacent layers allows for little error in FE predictions, (e.g., the 

modeled growth of a crack in one ply may either unduly arrest or promote crack growth in an 

adjacent ply if not well-predicted). More study is required to determine the significance of the 

link-up, which is likely to be specimen and loading dependent. It is anticipated that a predicted 
36 

delay in link-up will be significant for fatigue loading conditions. Analyses are currently in 

progress to assess the capabilities of refinements to the in-house continuum damage mechanics 

(CDM) finite element (FE) subroutine in properly representing these interactions. 

 45/90 interface 90/-45 interface -45/0 interface 

(a) CDM Predictions 

 45/90 interface 90/-45 interface -45/0 interface 

(b) X-FEM Predictions 

 45/90 interface 90/-45 interface -45/0 interface 

(c) X-Ray/CT Data 

Fig. 24. Comparison of analysis predictions using CDM [59] and X-FEM [14] models 

with X-Ray/CT data at a point just beyond peak load for [454/904/-454/04]s

specimens with a 1/8-inch diameter hole. 

Finally, experimental methods and data reduction methods are being developed to quantify 

crack densities and lengths at various load levels, to obtain crack initiation data, and to obtain 

crack-length versus load data. Crack initiation data for subsurface plies is being obtained, within 

the limits of current resolution, using acoustic emission coupled with X-Ray/CT. Crack initiation 

and crack growth data for surface plies are being obtained using digital image correlation. These 
37 

data will be used to aid in the calibration of the cohesive laws used in the analysis predictions, 

and to aid in the critical assessment of the analysis models. 

5 Summary 

Damage models for composites have quickly been increasing in diversity and complexity, and 

continuous technological breakthroughs are improving their predictive capabilities. Essential 

contributions to the development of numerous new modeling capabilities and to the formulation 

of the conditions under which given models can be expected to work have recently been made. 

An overview of recent developments in damage modeling for laminated composites was 

presented with the goal of defining the state of the art, the immediate technological barriers and 

for outlining potential solutions for research. In particular, the capabilities of advanced continuum 

damage mechanics models, cohesive zone models, and X-FEM models were reviewed and their 

pathological deficiencies were discussed. Issues of objectivity of the fracture propagation with the 

various models were discussed and the application of the methods was illustrated with examples 

of structural analysis. Finally, advances in experimental methods that are necessary for 

understanding the damage mechanisms in composites and for developing and validating the 

models were discussed. Non-destructive evaluation methods such as non-immersion UT scanners, 

digital image correlation, and X-Ray CT are providing rich experimental information on the 

response and failure of composite materials that is necessary for the development and validation 

of the computational models. 

References 

1. Wright P, Moffat A, Sinclair I, Spearing SM (2010) High Resolution Tomographic Imaging 

and Modelling of Notch Tip Damage in a Laminated Composite. Compos Sci Technol. 

70.(10): 1444-1452. 

2. Mollenhauer D, Ward L, Iarve E, et al. (2012) Simulation of Discrete Damage in 

Composite Overheight Compact Tension Specimens. Comp A. 43.(10): 1667-1679. 

3. van der Meer FP, Oliver C, Sluys LJ (2010) Computational Analysis of Progressive Failure 

in a Notched Laminate Including Shear Nonlinearity and Fiber Failure. Compos Sci 

Technol. 70.(4): 692-700. 

4. Fang XJ, Zhou ZQ, Cox BN, Yang QD (2011) High-Fidelity Simulations of Multiple 

Fracture Processes in a Laminated Composite in Tension. J Mech Phys Solids. 59.(7): 

1355-1373. 

5. González C, Llorca J (2007) Mechanical Behavior of Unidirectional Fiber-Reinforced 

Polymers under Transverse Compression: Microscopic Mechanisms and Modeling. 

Compos Sci Technol. 67.(13): 2795-2806. 

6. Maimí P, Camanho PP, Mayugo JA, Dávila CG (2007) A Continuum Damage Model for 

Composite Laminates: Part I - Constitutive Model. Mech Mater. 39.(10): 897-908. 

7. Camanho PP, Dávila CG, de Moura MFSF (2003) Numerical Simulation of Mixed-Mode 

Progressive Delamination in Composite Materials. J Compos Mater. 37.(16): 1415-1438. 

8. Krueger R (2002) The Virtual Crack Closure Technique: History, Approach and 

Applications. NASA/CR-2002-211628. , Hampton, VA. 

9. Mabson G, Lyle RD, Dopker B, Hoyt DM, Baylor JS, Graesser DL (2007) Fracture 

Interface Elements for Static and Fatigue Analysis. Proc 16th International Conference on 

Composite Materials, Kyoto, Japan. 

10. van der Meer FP, Sluys LJ (2009) Continuum Models for the Analysis of Progressive 

Failure in Composite Laminates. J Compos Mater. 43.(20): 2131-2156. 
38 

11. Iarve EV, Mollenhauer D, Kim R (2005) Theoretical and Experimental Investigation of 

Stress Redistribution in Open Hole Composite Laminates Due to Damage Accumulation. 

Comp A. 36.(2): 163-171. 

12. CMH-17 (2011) Materials Usage, Design and Analysis, Vol. 2, Rev G. Composite 

Materials Handbook. SAE, Philadelphia, PA. 

13. van der Meer FP, Sluys B, Hallett SR, Wisnom M (2012) Computational Modeling of 

Complex Failure Mechanisms in Laminates. J Compos Mater. 46.(5): 603-623. 

14. Swindeman MJ, Iarve EV, Brockman RA, Mollenhauer DH, Hallett SR (2011) Strength 

Prediction in Open Hole Composite Laminates by Using Discrete Damage Modeling. 52nd 

AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, 

Denver, CO. 

15. Griffith AA (1921) The Phenomena of Rupture and Flow in Solids. Phil Trans Roy Soc of 

London. A 221. 163-197. 

16. Dávila CG, Rose CA, Camanho PP (2009) A Procedure for Superposing Linear Cohesive 

Laws to Represent Multiple Damage Mechanisms in the Fracture of Composites. Int J 

Fracture. 158.(2): 211-223. 

17. Spearing SM, Beaumont PWR, Ashby MF (1992) Fatigue Damage Mechanics of 

Composite Materials. II: A Damage Growth Model. Compos Sci Technol. 44.(2): 169-177. 

18. ASTM D5528-01 (2002) Standard Test Method for Mode I Interlaminar Fracture 

Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites. Annual Book 

of Astm Standards. American Society for Testing and Materials, West Conshohocken, PA. 

19. Tamuzs V, Tarasovs S, Vilks U (2001) Progressive Delamination and Fiber Bridging 

Modeling in Double Cantilever Beam Composite Specimens. Eng Fract Mech. 68.(5): 513-

525. 

20. Gutkin R, Laffan ML, Pinho ST, Robinson P, Curtis PT (2011) Modelling the R-Curve 

Effect and Its Specimen-Dependence. Int J Solids Structures. 48.(11-12): 1767-1777. 

21. Sorensen L, Botsis J, Gmür T, Humbert L (2008) Bridging Tractions in Mode I 

Delamination: Measurements and Simulations. Compos Sci Technol. 68.(12): 2350-2358. 

22. Hansen AL, Lund E, Sørensen BF (2007) Simulation of Delamination Including Fiber 

Bridging Using Cohesive Zone Models. Proc ECCOMAS Thematic Conference on 

Mechanical Response of Composites P.P. Camanho et al., Porto, Portugal. 

23. Airoldi A, Dávila CG (2012) Identification of Material Parameters for Modelling 

Delamination in the Presence of Fibre Bridging. Compos Struct. 94.(11): 3240-3249. 

24. Högberg JL, Sørensen BF, Stigh U (2007) Constitutive Behaviour of Mixed Mode Loaded 

Adhesive Layer. Int J Solids Structures. 44.(25-26): 8335-8354. 

25. Turon A, Camanho PP, Costa J, Renart J (2010) Accurate Simulation of Delamination 

Growth under Mixed-Mode Loading Using Cohesive Elements: Definition of Interlaminar 

Strengths and Elastic Stiffness. Compos Struct. 92.(8): 1857-1864. 

26. Sarrado C, Turon A, Renart J, Urresti I (2012) Assessment of Energy Dissipation During 

Mixed-Mode Delamination Growth Using Cohesive Zone Models. Comp A. 43.(11): 2128-

2136. 

27. Leone AF, Girolamo D, Dávila CG (2012) Progressive Damage Analysis of Bonded 

Composite Joints. NASA/TM–2012-217790, Hampton. 

28. Turon A, Costa J, Camanho PP, Dávila CG (2007) Simulation of Delamination in 

Composites under High-Cycle Fatigue. Comp A. 38.(11): 2270-2282. 

29. Turon A, Camanho PP, Costa J, Dávila CG (2006) A Damage Model for the Simulation of 

Delamination in Advanced Composites under Variable-Mode Loading. Mech Mater. 

38.(11): 1072-1089. 
39 

30. Dávila CG, Camanho PP, Turon A (2008) Effective Simulation of Delamination in 

Aeronautical Structures Using Shells and Cohesive Elements. J Aircraft. 45.(2): 663-672. 

31. Ladevèze P, LeDantec E (1992) Damage Modelling of the Elementary Ply for Laminated 

Composites. Compos Sci Technol. 43.(3): 257-267. 

32. Matzenmiller A, Lubliner J, Taylor RL (1995) A Constitutive Model for Anisotropic 

Damage in Fiber-Composites. Mech Mater. 20.(2): 125-152. 

33. Kachanov LM (1986) Introduction to Continuum Damage Mechanics. Martinus Nijhoff, 

Dordrech, NL. 

34. Lemaitre J, Desmorat R, Sauzay M (2000) Anisotropic Damage Law of Evolution. Euro J 

Mech A. 19.(2): 187-208. 

35. Dassault Systèmes (2011) Abaqus® 6.11 Documentation. Simulia Corp., Providence, RI, 

USA. 

36. Maimí P, Camanho PP, Mayugo JA, Dávila CG (2006) A Thermodynamically Consistent 

Damage Model for Advanced Composites. NASA-TM-2006-214282, Hampton, VA. 

37. Soden PD, Hinton MJ, Kaddour AS (1998) A Comparison of the Predictive Capabilities of 

Current Failure Theories for Composite Laminates. Compos Sci Technol. 58.(7): 1225-

1254. 

38. Hinton MJ, Kaddour AS, Soden PD (2002) A Comparison of the Predictive Capabilities of 

Current Failure Theories for Composite Laminates, Judged against Experimental Evidence. 

Compos Sci Technol. 62.(12-13): 1725-1797. 

39. Dávila CG, Camanho PP, Rose CA (2005) Failure Criteria for Frp Laminates. J Compos 

Mater. 39.(4): 323-345. 

40. Pinho ST, Dávila CG, Camanho PP, Iannucci L, Robinson P (2005) Failure Models and 

Criteria for Frp under in-Plane or Three-Dimensional Stress States Including Shear NonLinearity. NASA/TM-2005-213530, Hampton, VA. 

41. de Borst R (2008) Challenges in Computational Materials Science: Multiple Scales, MultiPhysics and Evolving Discontinuities. Comput Mater Sci. 43.(1): 1-15. 

42. Bažant ZP, Oh B (1983) Crack Band Theory for Fracture of Concrete. Mater Struct. 16.(3): 

155-177. 

43. Maimí P, Camanho PP, Mayugo JA, Dávila CG (2007) A Continuum Damage Model for 

Composite Laminates: Part II - Computational Implementation. Mech Mater. 39.(10): 909-

919. 

44. Hansen P, Martin R (1999) DCB, 4ENF and MMB Delamination Characterisation of 

S2/8552 and IM7/8552. N68171-98-M-5177. United States Army European Research 

Office of the US Army, London, England. 

45. Czabaj MW, Ratcliffe JG (2012) Comparison of Intralaminar and Interlaminar Mode-I 

Fracture Toughness of Unidirectional IM7/8552 Graphite/Epoxy Composite. Proc 

American Society for Composites 27th Technical Conference, Arlington, TX. 

46. Pinho ST, Robinson P, Iannucci L (2006) Fracture Toughness of the Tensile and 

Compressive Fibre Failure Modes in Laminated Composites. Compos Sci Technol. 66.(13): 

2069-2079. 

47. Catalanotti G, Camanho PP, Xavier JC, Dávila CG, Marques AT (2010) Measurement of 

Resistance Curves in the Longitudinal Failure of Composites Using Digital Image 

Correlation. Compos Sci Technol. 70.(13): 1986-1993. 

48. Parvizi A, Garrett K, Bailey J (1978) Constrained Cracking in Glass Fibre-Reinforced 

Epoxy Cross-Ply Laminates. J Mat Sci. 13. 195-201. 

49. Wang ASD (1984) Fracture Mechanics of Sublaminate Cracks in Composite Materials. 

Comp Tech Rev. 6.(2): 45-62. 
40 

50. Camanho PP, Dávila CG, Pinho ST, Iannucci L, Robinson P (2006) Prediction of in Situ 

Strengths and Matrix Cracking in Composites under Transverse Tension and in-Plane 

Shear. Comp A. 37.(2): 165-176. 

51. Hahn HT, Tsai SW (1973) Nonlinear Elastic Behavior of Unidirectional Composite 

Laminae. J Compos Mater. 7.(1): 102-118. 

52. Schuecker C, Dávila CG, Pettermann HE (2008) Modeling the Non-Linear Response of 

Fiber-Reinforced Laminates Using a Combined Damage/Plasticity Model. NASA/TM-2008-

215314, Hampton, VA. 

53. Ramberg W, Osgood WR (1943) Description of Stress-Strain Curves by Three Parameters. 

NACA-TN-902, Hampton. 

54. Sihn S, Kim RY, Kawabe K, Tsai SW (2007) Experimental Studies of Thin-Ply Laminated 

Composites. Compos Sci Technol. 67.(6): 996-1008. 

55. Green BG, Wisnom MR, Hallett SR (2007) An Experimental Investigation into the Tensile 

Strength Scaling of Notched Composites. Comp A. 38.(3): 867-878. 

56. Jirásek M (1998) Nonlocal Models for Damage and Fracture: Comparison of Approaches. 

Int J Solids Structures. 35.(31-32): 4133-4145. 

57. Laš V, Zemčík R (2008) Progressive Damage of Unidirectional Composite Panels. J 

Compos Mater. 42.(1): 25-44. 

58. Song K, Li Y, Rose CA (2011) Continuum Damage Mechanics Models for the Analysis of 

Progressive Failure in Open-Hole Tension Laminates. 52nd AIAA/ASME/ASCE/AHS/ASC 

Structures, Structural Dynamics and Materials Conference, Denver, CO. 

59. Leone F, Dávila CG (2011) Application of Mixed-Mode Cohesive Damage Laws to 

Continuum Damage Mechanics. In: Ferreira A (ed) Proc ICCS16 -16th International 

Conference on Composite Structures, Porto, Portugal. 

60. Camanho PP, Bessa MA, Catalanotti G, Vogler M, Rolfes R (2012) Modeling the Inelastic 

Deformation and Fracture of Polymer Composites - Part II: Smeared Crack Model. 

Mechanics of Materials. 

61. Moës N, Dolbow J, Belytschko T (1999) A Finite Element Method for Crack Growth 

without Remeshing. Int J Num Meth Eng. 46.(1): 131-150. 

62. Belytschko T, Parimi C, Moës N, Sukumar N, Usui S (2003) Structured Extended Finite 

Element Methods for Solids Defined by Implicit Surfaces. Int J Num Meth Eng. 56.(4): 

609-635. 

63. Rudraraju S, Salvi A, Garikipati K, Waas AM (2012) Predictions of Crack Propagation 

Using a Variational Multiscale Approach and Its Application to Fracture in Laminated 

Fiber Reinforced Composites. Compos Struct. 94.(11): 3336-3346. 

64. Huynh DBP, Belytschko T (2009) The Extended Finite Element Method for Fracture in 

Composite Materials. Int J Num Meth Eng. 77.(2): 214-239. 

65. van der Meer F, Sluys L (2009) A Phantom Node Formulation with Mixed Mode Cohesive 

Law for Splitting in Laminates. Int J Fracture. 158.(2): 107-124. 

66. Ling D, Yang Q, Cox B (2009) An Augmented Finite Element Method for Modeling 

Arbitrary Discontinuities in Composite Materials. Int J Fracture. 156.(1): 53-73. 

67. Iarve EV (2003) Mesh Independent Modelling of Cracks by Using Higher Order Shape 

Functions. Int J Num Meth Eng. 56.(6): 869-882. 

68. Patzák B, Jirásek M (2003) Process Zone Resolution by Extended Finite Elements. Eng 

Fract Mech. 70.(7-8): 957-977. 

69. Benvenuti E (2008) A Regularized Xfem Framework for Embedded Cohesive Interfaces. 

Comp meth Appl Mech Eng. 197.(49-50): 4367-4378. 
41 

70. Oliver J, Huespe AE, Sánchez PJ (2006) A Comparative Study on Finite Elements for 

Capturing Strong Discontinuities: E-Fem Vs X-Fem. Comp meth Appl Mech Eng. 195.(37-

40): 4732-4752. 

71. Iarve EV, Gurvich MR, Mollenhauer DH, Rose CA, Dávila CG (2011) Mesh-Independent 

Matrix Cracking and Delamination Modeling in Laminated Composites. Int J Num Meth 

Eng. 88.(8): 749-773. 

72. Swindeman MJ (2011), "A Regularized Extended Finite Element Method for Modeling the 

Coupled Cracking and Delamination of Composite Materials," Doctor of Philosophy in 

Mechanical Engineering, The School of Engineering, University of Dayton, Dayton, OH. 

73. Li X, Hallett SR, Wisnom MR, Zobeiry N, Vaziri R, Poursartip A (2009) Experimental 

Study of Damage Propagation in over-Height Compact Tension Tests. Comp A. 40.(12): 

1891-1899. 

74. Fang XJ, Yang QD, Cox BN, Zhou ZQ (2011) An Augmented Cohesive Zone Element for 

Arbitrary Crack Coalescence and Bifurcation in Heterogeneous Materials. Int J Num Meth 

Eng. 88.(9): 841-861. 

75. van der Meer FP, Sluys LJ (2010) Mesh-Independent Modeling of Both Distributed and 

Discrete Matrix Cracking in Interaction with Delamination in Composites. Eng Fract 

Mech. 77.(4): 719-735. 

76. van der Meer FP, Dávila CG (2013) Cohesive Modeling of Transverse Cracking in 

Laminates under in-Plane Loading with a Single Layer of Elements Per Ply. Int J Solids 

Structures. 50.(20–21): 3308-3318. 

77. Seebo J, Johnston P (2010) High Speed Position Synchronized Data Acquisition System for 

in-Situ Ultrasonic Inspection. Proc 2010 ASNT Spring Conference, Williamsburg, VA. 

78. Johnston PH, Appleget CD, Odarczenko MT (2013) Characterization of Delaminations and 

Transverse Matrix Cracks in Composite Laminates Using Multiple-Angle Ultrasonic 

Inspection. Proc AIP Conference Proceedings, 1011. 

79. Correlated Solutions (2012) VIC-3D User Manual, Columbia, SC. 

80. Leone FA, Girolamo D, Dávila CG (2012) Progressive Damage Analysis of Bonded 

Composite Joints. NASA/TM-2012-217790, Hampton, VA. 

81. Scott AE, Mavrogordato M, Wright P, Sinclair I, Spearing SM (2011) In Situ Fibre 

Fracture Measurement in Carbon-Epoxy Laminates Using High Resolution Computed 

Tomography. Compos Sci Technol. 71.(12): 1471-1477. 

82. Hallett SR, Green BG, Jiang WG, Wisnom MR (2009) An Experimental and Numerical 

Investigation into the Damage Mechanisms in Notched Composites. Comp A. 40.(5): 613-

624. 
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

Analysis Methods for Progressive Damage of Composite Structures

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Rose, Cheryl A.; Davila, Carlos G.; Leone, Frank A.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20257

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 24

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

This document provides an overview of recent accomplishments and lessons learned in the development of general progressive 

damage analysis methods for predicting the residual strength and life of composite structures. These developments are 

described within their State-of-the-Art (SoA) context and the associated technology barriers. The emphasis of the authors is on

developing these analysis tools for application at the structural level. Hence, modeling of damage progression is undertaken at

the mesoscale, where the plies of a laminate are represented as a homogenous orthotropic continuum. The aim of the present 

effort is establish the ranges of validity of available models, to identify technology barriers, and to establish the foundations of 

the future investigation efforts. Such are the necessary steps towards accurate and robust simulations that can replace some of

the expensive and time-consuming “building block” tests that are currently required for the design and certification of 

aerospace structures.

15. SUBJECT TERMS

Composites; Damage models; Material characterization; Nondestructive evaluation

18. NUMBER

 OF 

 PAGES

48

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

794072.02.07.03.03

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TM-2013-218024

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 07 - 2013