Evan J. Pineda

Glenn Research Center, Cleveland, Ohio

Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

Numerical Implementation of a Multiple-ISV 

Thermodynamically-Based Work Potential 

Theory for Modeling Progressive Damage 

and Failure in Fiber-Reinforced Laminates

NASA/TM—2011-217401

December 2011
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 Reports 

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 peer-reviewed formal 

professional papers but has 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. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

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 via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
Evan J. Pineda

Glenn Research Center, Cleveland, Ohio

Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

Numerical Implementation of a Multiple-ISV 

Thermodynamically-Based Work Potential 

Theory for Modeling Progressive Damage 

and Failure in Fiber-Reinforced Laminates

NASA/TM—2011-217401

December 2011

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

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

Trade names and trademarks are used in this report for identification 

only. Their usage does not constitute an official endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

Level of Review: This material has been technically reviewed by technical management. 
1 Introduction

Predictive capabilities of numerical tools for progressive damage and failure analysis (PDFA) of composite structures are predicated upon the robustness, accuracy,

and objectivity of the tools. As sophisticated numerical tools become more feasible

for widespread use in industry, significant weight and cost saving in the design of

lightweight composite structures will be apparent. Virtual testing of materials, incorporating PDFA, can be used to evaluate the viability of materials or configurations

prior to further scrutiny via physical testing. This grants analysts more flexibility during preliminary structural design stages and will ultimately manifest as more efficient

and cost effective designs.

Continuum damage mechanics (CDM) has emerged as a viable option for predicting the non-linear behavior of composite structures. The first CDM theory was developed by Kachanov (1958, 1986). Subsequently, many publications on this subject

were produced, including numerous books [Talreja (1985a); Lemaitre and Chaboche

(1994); Lemaitre (1996); Krajcinovic (1996); Voyiadjis and Kattan (2005)]. Typically, a set of scalar damage variables, or internal state variables (ISVs), introduce

anisotropic damage into the composite constituent behavior by penalizing the components of the material stiffness tensor, and non-linear functions are used to control

the damage evolution. Various authors have used crack density, geometry, strain

energy release rate, and other crack features to characterize the damage evolution

[Dvorak et al. (1985); Talreja (1985b); Laws and Dvorak (1988); Lee et al. (1989);

Nairn (1989); Tan and Nuismer (1989); Gudmundson and Ostlund  ̈ (1992); McCartney (1992a, 1998)]. Others postulate damage evolution laws and characterize those

laws using experiments [Allen et al. (1987a,b); Talreja (1994); Paas et al. (1992);

Matzenmiller et al. (1995); Bednarcyk et al. (2010)]. CDM models must also employ

failure criteria to indicate damage initiation. More recently, increasingly sophisticated failure criteria have been developed to better represent the phenomenological

NASA/TM—2011-217401 1
behavior of a damaging composite lamina [Puck and Sch ̈urmann (1998, 2002); Pinho

et al. (2005)] and used in conjunction with CDM.

CDM techniques offer computationally efficient and readily implementable means

to capturing the effects of damage and failure in composite materials. Unfortunately,

the majority of the criteria and evolution laws are formulated upon phenomenological

observations of the directional dependence of damage evolution, rather than modeling

the physics of the actual damage mechanisms. Separate damage variables are used to

degrade different components of the stiffness tensor (directions) depending on whether

the damage is said to accrue in the matrix or fiber constituents of the composite, but

the variables do not explicitly distinguish between the separate damage mechanisms.

Furthermore, many theories involve a multitude of parameters that are difficult to

measure and must be calibrated to correlate with experimental data.

When implemented within the finite element method (FEM), many PDFA methodologies that utilize CDM breakdown when the material enters the post-peak strain

softening regime locally within an element. Loss of positive definiteness of the tangent

stiffness tensor leads to pathological mesh dependence [Baˇzant and Cedolin (1979);

Pietruszczak and Mroz (1981)]. To overcome this deficiency, Baˇzant (1982) developed the smeared crack, or crack band, model that introduces a characteristic element length into the formulation of the damage evolution. The original formulation

assumed that the mode I crack band always aligns with the principle axes [Baˇzant

and Oh (1983)]. de Borst and Nauta (1985) altered the formulation to accommodate a fixed crack band under mixed mode conditions. An encompassing overview of

smeared crack band models is provided by Spencer (2002).

A thermodynamically-based, work potential theory, known as Schapery theory

(ST), was developed for modeling matrix microdamage in fiber-reinforced laminates

(FRLs) [Lamborn and Schapery (1988); Schapery (1989, 1990); Lamborn and Schapery

(1993)]. Sicking (1992) and Schapery and Sicking (1995) extended the formulation

NASA/TM—2011-217401 2
to include the effects of transverse cracking by adding an additional internal state

variable (ISV) and predicted the evolution of microdamage and transverse cracking

in coupon laminates analytically. Pineda et al. (2010) implemented this extended

formulation in a numerical setting to simulate the failure of a buffer strip-reinforced,

center-notched panel (CNP). However, due to the cumbersome nature of the evolution equations, the microdamage and transverse cracking evolution equations were

decoupled to arrive at a more efficient implementation. Since no characteristic length

is introduced into the formulation, the theory produces mesh-dependent results in a

computational setting.

The ST formulation is modified here, resulting in the enhanced Schapery theory

(EST), to include the effects of macroscopic transverse and shear matrix cracking, as

well as fiber breakage, using an approach that differs from Sicking (1992); Schapery

and Sicking (1995); Pineda et al. (2010). A deliberate distinction between damage and

failure is made. Damage is defined as the effects of any structural changes resulting in

a non-linear response that preserves the positive definiteness of the tangent stiffness

tensor of the material. Conversely, failure is considered to be the consequence of

structural changes that cause post-peak strain softening in the stress versus strain

response of the material. Here, matrix microdamage is categorized as a damage

mechanism, but macroscopic matrix cracking and fiber breakage are hypothesized to

be failure mechanisms resulting from damage localization. The traditional ISV used in

ST is maintained to model microdamage. Upon failure initiation, the element domain

is no longer considered a continuum, and a smeared crack approach is used to model

the embedded discontinuities. Three new ISVs, which incorporate the characteristic

length of the finite element, dictate the evolution of the failure mechanisms. The EST

formulation presented in Section 2 offers non-linear progressive damage coupled with

mesh objective, post-peak strain softening.

Mesh objectivity is demonstrated in Section 3. In Section 4, EST is verified against

NASA/TM—2011-217401 3
experimental results for two center-notched panels (CNPs). Global load versus deflection data, local strain gage data, as well as observed failure mechanisms obtained

from experiments performed at the NASA Langley Research Center (LaRC) and exhibited in Bogert et al. (2006); Satyanarayana et al. (2007) are compared to numerical

results.

2 Enhanced Schapery Theory

The previously developed ST ([Schapery (1990, 1995); Schapery and Sicking (1995);

Basu et al. (2006); Pineda et al. (2009, 2010)]) is extended to accommodate mesh

objective, post peak strain softening. Separate ISVs are used to govern the evolution

of matrix microdamage, transverse (mode I) matrix failure, shear (mode II) matrix

failure, and fiber breakage (mode I). The first and second laws of thermodynamics

are enforced, establishing thermodynamically consistent evolution laws for progressive

matrix microdamage, as well as post-peak failure. The following sections detail the

formulation of this work potential theory.

2.1 Thermodynamically-based Work Potential Framework

As a material is loaded, a measure of the work potential facilitates modeling structural

changes in the material, such as microcracking, which affect the elastic properties of

the material. Energy that is not dissipated is recovered when the structure is unloaded, and the magnitude of energy recovered is contingent upon the degraded,

elastic properties at the previously attained maximum strain state. It is assumed,

upon subsequent reloading, that the material behaves linearly, exhibiting the elastic properties observed during unloading, until the material reaches the preceding

maximum strain state. After this state is achieved, structural changes resume, affecting degradation of the instantaneous elastic moduli of the material. This process

NASA/TM—2011-217401 4
is shown in the uniaxial stress-strain curve displayed in Figure 1. The shaded area

above the unloading line represents total dissipated potential WS, and the triangular

area underneath is the total elastic strain energy density WE. It is assumed that the

material behaves as a secant material and there is no permanent deformation upon

unloading. This a reasonable assumption for FRLs Sicking (1992); however, plastic

deformation can also be incorporated, if necessary [Schapery (1990)]. Extension to

treat viscoelastic and viscoplastic response is outlined in Hinterhoelzl and Schapery

(2004).

Both WE and WS are functions of a set of ISVs, Sm,(m = 1, 2, M). These ISVs

account for any inelastic structural changes in the material. Differentiating WS with

respect to any ISV Sm, assuming limited path-dependence [Schapery (1990)], yields

the thermodynamic force, fm, available for advancing structural changes associated

with the mth ISV.

fm = ∂Ws

∂Sm

(1)

It is shown in Schapery (1989, 1990) that the total work potential is stationary with

respect to each ISV.

∂WT

∂Sm

= 0 (2)

Additionally, Rice (1971) utilized the second law of thermodynamics to establish the

inequality:

fmS ̇

m ≥ 0 (3)

which suggests that “healing” is not allowed for a material undergoing structural

changes. Equations (1), (2), and (3) form the foundation of a thermodynamicallybased work potential theory for modeling non-linear structural changes in a material

exhibiting limited path-dependence.

NASA/TM—2011-217401 5
2.2 Multiple ISV Formulation of ST to Account for Multiple

Damage and Failure Mechanisms

Due to the generality of the evolution equations, Equations (2) and (3), the work

potential theory can account for any number and type of structural changes that

may occur in a material. This is especially useful for modeling progressive damage

in composites because the heterogeneity of the composite, and multiaxiality of the

local fields, enables multiple damage mechanisms to arise during a typical loading

history. For instance in the matrix phase alone, microdamage accrues until its effects

are superseded by the growth of larger transverse cracks. Microdamage is considered

the advancement of microcracks, voids, fissures, shear bands, and other flaws that are

present in the matrix of a composite [Sicking (1992); Schapery and Sicking (1995);

Basu et al. (2006); Ng et al. (2010)]. The size of these flaws is typically on the order

of that of the fiber or smaller. Transverse cracks nucleate from pre-existing flaws

within the matrix but grow parallel to the fibers and span the thickness of the lamina

[Talreja (1985b); Allen et al. (1987a); Laws and Dvorak (1988); Gudmundson and

Ostlund  ̈ (1992); Yang and Cox (2005); Noda et al. (2006); Green et al. (2007)]. Often,

the growth of individual transverse cracks is extremely rapid; however, the effects of

transverse cracking on the stiffness of a composite laminate can be progressive if

multiple cracks form over an extended period of time and throughout an expansive

volume. Eventually, transverse cracking is succeeded by more catastrophic damage

mechanisms including interlaminar delamination, fiber breakage, pullout and bridging

associated with macroscopic laminate fracture [McCartney (1992a,b); Hallett et al.

(2008)].

The present EST formulation assumes that three major intralaminar mechanisms

are responsible for all observed non-linearities in the stress-strain curve of a composite

lamina: matrix microdamage, matrix macroscopic cracking, and axial fiber failure.

Each of these mechanisms can be accommodated by partitioning the total dissipated

NASA/TM—2011-217401 6
energy density, WS, into portions associated with each mechanism.

Matrix microdamage is the primary cause of observed non-linearity in the stress

versus strain response of many polymer matrix composites (PMCs) (i.e. systems exhibiting negligible non-linear elasticity, plasticity or viscous effects) up to localization

of microdamage into more severe failure mechanisms, such as transverse cracking,

fiber breakage, kink band formation, or delamination. Microdamage can be considered the combination of matrix microcracking, micro-void growth, shear banding,

and fiber-matrix debonding. Figure 2 shows a typical uniaxial response of a material

exhibiting microdamage evolution, where the recoverable energy potential is given

by W and the potential dissipated into evolving structural changes associated with

microdamage is given by S.

Typically, matrix microdamage continues to grow until the onset of more catastrophic failure mechanisms initiate. It should be noted that this work explicitly

distinguishes between damage and failure in the following manner:

Damage - Structural changes in a material that manifest as pre-peak non-linearity in

the stress-strain response of the material through the degradation of the secant

moduli.

Failure - Structural changes that result from damage localization in a material and

manifest as post-peak strain softening in the stress-strain response of the material.

Here, three major failure mechanisms, which are distinct from the microdamage mode,

are considered: transverse (mode I) matrix cracking, shear (mode II) matrix cracking,

and axial (mode I) fiber fracture. These failure modes are consistent with the in-plane

failure typically observed in PMC laminates. It is assumed that the evolution of these

mechanisms yields an immediate reduction in the load-carrying capability of a local

subvolume where the mechanism is active. Three ISVs are used to account for mode

NASA/TM—2011-217401 7
I matrix cracking, mode II matrix cracking, and mode I fiber failure, respectively:

Sm

I , SmII , and Sf

I . These ISVs are defined completely in Section 2.3, and are taken to

be the potentials required to advance structural changes associated with these failure

mechanisms.

At any given state the total dissipated energy density WS can be calculated as a

sum of energy dissipated through the aforementioned damage and failure mechanisms,

given by the four ISVs.

WS = S + Sm

IF + SmIIF + Sf

I (4)

According to the first law of thermodynamics, the total work potential (ignoring

thermal dissipation) is given by the sum of the elastic strain energy density and the

potentials associated with each of the damage or failure mechanisms.

WT = WE + S + Sm

I + SmII + Sf

I (5)

where WE is the elastic strain energy density. Invoking the stationarity principle,

Equation (2),

∂WE

∂S = −1

∂WE

∂Sm

I

= −1

∂WE

∂Sm

II

= −1

∂WE

∂Sf

I

= −1

(6)

NASA/TM—2011-217401 8
and the Second Law of Thermodynamics, Equation (3), gives:

S ̇ ≥ 0

S ̇ m

IF ≥ 0

S ̇ m

IIF ≥ 0

S ̇ f

F ≥ 0

(7)

Equations (6) and (7) constitute the evolution equations for damage and failure in

a material associated with matrix microdamage, matrix cracking, and fiber breakage

in tension.

It should be noted, that EST can also account for kink band formation under

axial compression [Schapery (1995); Basu (2005); Basu et al. (2006)]; although, the

applied loading in the examples presented in Sections 3 and 4 are tensile, and kink

banding does not occur. As the lamina is loaded, the fibers in the composite rotate

by some angle φ, given by the deformation gradient in the model. To model the kink

band mechanism, all calculations are then executed in the instantaneous fiber frame

given by φ; therefore, fibers rotation induces larger shear strains, γ12. Increased shear

strain yields more damage, leading to a reduction in the shear modulus. The increase

in shear compliance allows for further progression of the shear strain. Under axial

compression, this leads to a runaway instability, and a kink band will form.

2.3 Failure Initiation

Matrix microdamage requires no initiation criterion. For low strain levels, the microdamage ISV S remains small and its effects on the composite moduli are not apparent.

As S evolves, with increased strains, its effects on the stress-strain response of the

NASA/TM—2011-217401 9
composite become more noticeable. However, it is postulated that the evolution of the

failure mechanisms immediately yield a negative tangent stiffness; therefore, initiation criteria are required. Furthermore, criteria are required to mark failure initiation

because the macroscopic cracks responsible for failure may result from localization

of microdamage, or they may nucleate from pre-existing flaws in the material not

necessarily associated with microdamage.

EST is implemented in homogenized laminae; therefore, phenomenological criteria

must be utilized that account for the composite microstructure. The Hashin-Rotem

failure criterion incorporates separate equations for matrix failure and fiber failure

initiation [Hashin and Rotem (1973)]. The matrix failure criterion involves contributions from both the transverse (22) and shear (γ12) strains.

22

YT

2

+

γ12

Z

2

= 1 22 ≥ 0

22

YC

2

+

γ12

Z

2

= 1 22 < 0

(8)

where YT is the transverse lamina failure strain in tension, YC is the transverse failure

lamina strain in compression, and Z is the shear failure strain. The fiber failure

criterion only involves the axial strain 11.

11

XT

2

= 1 11 ≥ 0 (9)

where XT is the maximium allowable axial strain of the lamina. A local, lamina

coordinate frame is chosen such that, 1- is the axial direction of the fibers, 2- is the

in-plane transverse direction, and 3- is the out-of-plane direction. When Equation

(8) is satisfied the matrix failure ISVs Sm

I and SmII are activated, and when Equation

(9) is satisfied fiber failure evolution Sf

I is permitted; otherwise, S remains the only

active ISV. Upon satisfaction of either Equation (8) or Equation (9), it is assumed

NASA/TM—2011-217401 10
that the more severe failure mechanisms dominate, superseding the effects of matrix

microdamage; therefore, S ̇ = 0, and additional microdamage is precluded.

2.4 Use of Traction-Separation Relationships to Define the

Failure Potentials

Sicking (1992) and Schapery and Sicking (1995) used a single ISV to model the effects

of transverse cracking on a composite lamina. Similar to microdamage, the transverse

and shear moduli were related to transverse crack evolution through a set damage

functions obtained from coupon experiments. Predictions of the non-linear response

of numerous laminates were presented assuming a homogenous strain state in the

laminates. Pineda et al. (2010) implemented the dual-ISV formulation of ST for

predicting microdamage and transverse cracking within FEM to model the response

of a center-notched laminate that was reinforced with buffer strips. The original

formulation required the solution of two, coupled, bi-variate polynomials, which in an

FEM framework became extremely computationally intensive. Thus, Pineda et al.

(2010) decoupled the microdamage and transverse cracking evolution equations.

In the aforementioned publications, it was assumed that the transverse cracking

affected the relationships between stress and strain. However, the existence of a

macroscopic crack invalidates the assumption of a continuum. Here, it is presumed

that failure arises from the evolution of cohesive cracks within the continuum, and the

ISVs associated with failure (axial, transverse, and shear) influence the relationship

between traction on the crack faces and the crack-tip opening displacement. The

satisfaction of Equations (8) and/or (9) indicates the material behavior transitions

from that of a damaging continuum to that of a cohesive crack, and the essential

fields become traction and separation, rather than stress and strain (see Figure 3).

Once a cohesive crack initiates in the continuum, opening of the crack yields a

reduction in traction on the crack faces at the crack tip. If subsequently the crack

NASA/TM—2011-217401 11
is closed, it is assumed that traction at the crack tip will unload linearly towards

the origin of the traction versus separation law (see Figure 4). The strain energy

release rate (SERR) Gj

M is taken as the total energy dissipated per unit area of new

surface that is created through crack advancement and can be calculated as the area

under the traction-separation law (for a given traction and separation pair) minus

the energy per area that can potentially be recovered by unloading.

Gj

M =

 δj

M

0

t

j

MdδjM − 1

2

t

j

Mδ

j

M (10)

where j indicates the material (fiber f or matrix m), M represents the corresponding

mode (mode I or mode II), δ

j

M is the crack tip opening displacement in mode M and

material j, and t

j

M is the corresponding traction at the crack tip.

Theoretically, the shape of the traction-separation laws for mode I crack growth

in the fiber, and mode I and II crack growth in the matrix can take any shape.

Gustafson and Waas (2009) investigated triangular, trapezoidal, beta distribution,

and sinusoidal traction-separation laws in discrete cohesive zone method (DCZM) elements and determined that the shape only affected convergence of the FEM solver,

but not the overall results. For simplicity, it is assumed here that all three types of

cracks obey triangular traction-separation laws, presented in Figure 4. The total area

under the traction-separation curves is controlled by the corresponding material fracture toughness in the appropriate mode, where Gf

IC is the mode I fracture toughness

of the fiber, Gm

IC is the mode I fracture toughness of the matrix, and GmIIC is the mode

II fracture toughness of the matrix. The cohesive strengths of the materials t

f

IC (mode

I fiber strength), t

m

IC (mode I matrix strength), and t

m

IIC (mode II matrix strength)

are given by the stresses in the continuum when Equations (8) and/or (9) are satisfied. Mode I, normal cracks are not allowed to grow under compression, but mode II,

shear cracks can evolve under normal compression. Therefore, the mode I tractionNASA/TM—2011-217401 12
separation laws for the fiber and matrix (Figures 4a and 4b) do not accommodate

negative crack tip displacements. However under negative mode II displacement (see

Figure 4c), the traction on the crack faces will increase linearly until the maximum,

previously attained displacement magnitude is reached, after which, the crack faces

will resume unloading according to the negative portion of the traction-separation

law. The traction-separation laws exhibited in Figure 4 do not require any initial,

fictitious, pre-peak stiffness because the cracks are embedded within a continuum.

This is an advantage over the use of DCZM elements which do require an initial stiffness because these interfacial elements do not actually represent physical material

within the model and must attempt to simulate initially perfect bonding between

adjacent material domains [Xie et al. (2006); Gustafson (2008)]. If set incorrectly,

these fictitious stiffnesses can cause numerical problems [Turon et al. (2006)].

Although no mode I crack can advance under compression, it is possible for postpeak softening to occur under compressive loading situations. For instance a kink

band could form under global axial compression, or the matrix could fail in local

shear due to internal friction (Mohr-Coulomb) in quasi-brittle materials under transverse compression [Hoek and Bieniawski (1965)]. Since these failure mechanisms

involve local shear at a the fiber/matrix scales which is typically below the operating

lamina/laminate scale, it appears that these mechanisms evolve under mode I compression. In a model containing homogenized laminae, there is no subscale shear to

drive these compressive failure modes. However, EST could be extended further to

incorporate these mechanisms through phenomenological accessions. The methods,

developed by Basu (2005); Basu et al. (2006) and described in Section 2.3, can be used

to track the instantaneous fiber angle, and a critical fiber angle can be assigned to

indicate the initiation of post-peak softening due to kink band formation. Similarly,

a matrix compression criterion, such as the one developed by Puck and Sch ̈urmann

(1998, 2002) could be used to signal the initiation of Mohr-Coulomb compressive

NASA/TM—2011-217401 13
failure. The traction-separation laws for mode I fiber compression and mode I transverse matrix compression could be adjusted to include the post-peak softening effects

of microbuckling and Mohr-Coulomb matrix failure. These postulated, compressive,

mode I traction-separation laws could account for energy released through these subscale failure modes in a homogenous model at the lamina/laminate scale. However,

the examples presented in this chapter are tension dominated, and extension of the

theory to accommodate apparent mode I compressive failure is left for future work.

Using the traction-separation laws in Figure 4, the SERR can be calculated with

Equation (10).

Gf

I = 1

2

t

f

ICδ

f

I (11)

Gm

I = 1

2

t

m

ICδmI (12)

Gm

II = 1

2

t

m

IICδmII (13)

It is assumed that the energy released due to cracking is smeared over the entire

element [Baˇzant (1982); Baˇzant and Oh (1983)]. Thus, the dissipation potentials in

an element resulting from macroscopic cracking are related to the SERRs using the

suitable element dimensions.

Sf

I = Gf

I

l

(θ+90◦) e

(14)

Sm

I = Gm

I

l

(θ) e

(15)

Sm

II = Gm

II

l

(θ) e

(16)

If there is a single integration point in the element, l

(θ+90◦) e is the length of a line

running perpendicular to fiber direction in the element that intersects two edges of

the element and the integration point, and l

(θ) e is the length of a line that is parallel

to the fiber direction in the element that intersects two edges of the element and

the integration point. If there is more than one integration point in the element,

NASA/TM—2011-217401 14
the element can be partitioned into a number of subvolumes equal to the number of

integration points, and the lengths l

(θ+90◦) e and l(θ) e are lengths of lines that intersect

the corresponding integration point as well as two element edges or integration point

subvolume boundaries. Incorporating a length scale into the ISVs results in mesh

objective, post-peak, softening. This is elaborated upon further in Section 4.

2.5 EST Evolution Equations for a Fiber-Reinforced Lamina

To arrive at the evolution equations for the four ISVs, the elastic strain energy density

must be defined for a material which may contain cohesive cracks. Therefore, the

elastic strain energy WE is comprised of a contribution from the continuum W and

any possible cohesive cracks Wj

M. The plane stress, elastic strain energy density in

the continuum is defined as

W = 1

2

(E11

2

11 + E22(S)

2

22 + G12(S)γ212) + Q121122 (17)

where stress in the laminae are related to strain assuming plane stress conditions.

σ11 = Q1111 + Q1222

σ22 = Q1211 + Q2222

τ12 = Q66γ12

(18)

where γ12 is the engineering shear strain and

Q11 = E11

1 − ν12ν21

Q22 = E22

1 − ν12ν21

Q12 = ν12Q22

Q66 = G12

ν21 = ν12E22

E11

(19)

NASA/TM—2011-217401 15
where E11 is the axial elastic modulus, E22 is the transverse elastic modulus, ν12 is

the Poisson’s ratio, ν21 is the transverse Poisson’s ratio, and G12 is the elastic shear

modulus. After assuming that the quantity ν12ν21 << 1, Equations (19) simplify,

Q11 = E11

Q22 = E22

Q12 = ν12Q22

Q66 = G12

(20)

Note that only the transverse and shear moduli (E22 and G12) are functions of S

since matrix microdamage only accrues in the matrix of the laminae. The Poisson’s

ratio is assumed to evolve such that the quantity Q12 = E22ν12 remains constant;

however, this restriction can be relaxed if deemed necessary. The degraded moduli

are related to the virgin moduli (E220 and G120) and the ISV through a set of microdamage functions (es(S) and gs(S)) that are obtained from three uniaxial coupon

tests [Schapery (1989); Sicking (1992); Schapery and Sicking (1995)].

E22 = E220es(S) (21)

G12 = G120gs(S) (22)

Degrading E22 and G12 exclusively is consistent with the intralaminar damage typically observed in PMC laminates.

The elastic strain energy density of the cohesive cracks are defined as the recoverable energy per unit crack surface area smeared over the entire element.

Wf

I = t

f

I δ

f

I

2l

(θ+90◦) e

(23)

W m

I = t

m

I δmI

2l

(θ) e

(24)

NASA/TM—2011-217401 16
W m

II = t

m

II δm

II

2l

(θ) e

(25)

The tractions in Equations (23)-(25) can be related to the secant stiffness’ in the

traction-separation laws kj

M.

t

f

I = kfI δ

f

I (26)

t

m

I = kmI δmI (27)

t

m

II = kmII δmII (28)

Hence, the total elastic strain energy density in the continuum is given by

WE =1

2



E11

2

11 + E22(S)

2

22 + G12(S)γ212

+ Q121122

+

kf

I (SfI )

2l

(θ+90◦) e

δ

f

I

2

+

km

I (SmI )

2l

(θ) e

δm

I

2 +

km

II (SmII )

2l

(θ) e

δm

II

2

(29)

Substituting Equation (29) into Equations (6) gives the ISV evolution equations.

1

2





2

22E220

des

dSr

+ γ2

12G120

dgs

dSr



= −3S2

r (30)

1

2l

(θ+90◦) e

dkf

I

dSf

I

δ

f

I

2

= −1 (31)

1

2l

(θ) e

dkm

I

dSm

I

δm

I

2 = −1 (32)

1

2l

(θ) e

dkm

II

dSm

II

δm

II

2 = −1 (33)

The use of a reduced ISV Sr = S 1

3 has been employed in Equation (30). Sicking

(1992) has shown that the use of this reduced ISV yields polynomial forms of the

microdamage functions in Equations (21) and (22). Using the chain rule and the fact

that

dSf

I

dδf

I

= t

f

IC

2l

(θ+90◦) e

(34)

NASA/TM—2011-217401 17
dSm

I

dδm

I

= t

m

IC

2l

(θ) e

(35)

dSm

II

dδm

II

= t

m

IIC

2l

(θ) e

(36)

by Equations (11)-(16), the cohesive secant stiffnesses are determined.

kf

I = −

 t

f

IC

δ

f

I

2 dδf

I (37)

km

I = −

 t

m

IC

δm

I

2 dδm

I (38)

km

II = −

 t

m

IIC

δm

II

2 dδm

II (39)

Evaluating the integrals in Equations (37)-(39), while enforcing kj

M = 0 when δ

j

M =

2Gj

MC

t

j

MC

results in expressions for kj

M in terms of δ

j

M.

kf

I = t

f

IC 	

1

δ

f

I

− t

f

IC

2Gf

IC 


(40)

km

I = t

m

IC 1

δm

I

− t

m

IC

2Gm

IC



(41)

km

II = t

m

IIC 1

δm

II

− t

m

IIC

2Gm

IIC



(42)

The thermodynamically consistent stiffnesses derived in Equations (40)-(42) can also

be derived directly from the traction-separation laws using geometry.

Finally, it is assumed that following failure initiation the strains are related to the

crack tip opening displacements by

l

(θ+90◦) e 11 = l(θ+90◦) e C

11 + δ

f

I (43)

l

(θ)

e 22 = l

(θ)

e 

C

22 + δmI (44)

NASA/TM—2011-217401 18
l

(θ)

e γ12 = l

(θ)

e γC12 + 2δmII (45)

where C

11, C22, and γC12 are the strains when Equations (8) and/or (9) are satisfied.

Equations (43)-(45) imply that the strain in the continuum remains at the values

obtained when failure initiates, and that any incremental change in the global strain

after failure initiation is used wholly to advance the crack tip opening displacement.

To account for changes in the continuum strain after failure initiates, it can be assumed that the stress state in the cracked body is homogenous and the tractions on

the crack tip faces are equal to the stresses in the continuum. Then, the strains in

Equation (29) can be formulated in terms of the cohesive secant stiffnesses and the

crack tip opening displacement. However, it is assumed that the evolution of strain

in the continuum is negligible once cohesive cracks form. Equations (43)-(45) can be

utilized in Equations (40)-(42) to obtain kj

M as functions of the global strain at an

integration point.

kf

I = t

f

IC 

1

l

(θ+90◦) e (11 − C

11)

− t

f

IC

2Gf

IC 

(46)

km

I = t

m

IC 

1

l

(θ) e (22 − C

22)

− t

m

IC

2Gm

IC 

(47)

km

II = t

m

IIC 

2

l

(θ) e (γ12 − γC

12)

− t

m

IIC

2Gm

IIC 

(48)

Once failure initiates, the effects of failure supersede the effects of microdamage

and evolution of S ceases. The cohesive stiffness in a cracked element is calculated

using Equations (46)-(48) for a given strain state; then, Equations (26)-(28) and

(43)-(45) are used to calculate the tractions on the crack tip faces and the crack tip

opening displacement. It is assumed that the stress state in the integration point

subvolume of the element is homogenous, and the tractions on the crack tip faces are

equal to the stresses in the element. Lastly, the axial, transverse, and shear moduli

NASA/TM—2011-217401 19
of the element can be calculated [Baˇzant and Oh (1983)]:

E11 =

⎧

⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎩

1

E110

− 11 − C

11

t

f

IC 

1 +

l

(θ+90◦) e t

f

IC

2Gf

IC



11 − 

C

11



⎫

⎪⎪⎪⎪⎬

⎪⎪⎪⎪⎭

−1

(49)

E22 =

⎧

⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎩

1

E∗

22

− 22 − C

22

t

m

IC 

1 +

l

(θ) e t

m

IC

2Gm

IC



22 − 

C

22



⎫

⎪⎪⎪⎪⎬

⎪⎪⎪⎪⎭

−1

(50)

G12 =

⎧

⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎩

1

G∗

12

− γ12 − γC

12

2t

m

IIC 

1 +

l

(θ) e t

m

IIC

4Gm

IIC



γ12 − γC

12



⎫

⎪⎪⎪⎪⎬

⎪⎪⎪⎪⎭

−1

(51)

where E∗

22 and G∗12 are the degraded transverse and shear moduli, due to microdamage, when Equation (8) is satisfied.

For visualization purposed in the FEM simulations, degradation parameters are

defined which relate the current, degraded stiffnesses to their original values upon

failure initiation.

Df

I = 1 − E11

E∗

11

(52)

Dm

I = 1 − E22

E∗

22

(53)

Dm

II = 1 − G12

G∗

12

(54)

The degradation parameter can have a minimum value of zero, which indicates that

no degradation has occurred, or a maximum value of one, signaling that the corresponding modulus has been completely diminished.

The negative tangent stiffness of the stress-strain curve necessary for post-peak

strain softening to occur imposes a restriction the maximum allowable element size,

NASA/TM—2011-217401 20
as shown by Baˇzant and Oh (1983).

l

(θ+90◦) e <

2Gf

ICE11

t

f

IC

2 (55)

l

(θ)

e < min 2Gm

ICE∗

22

tm

IC

2 ,

2Gm

IICG∗

12

tm

IIC

2



(56)

The analyst must be careful to ensure the dimensions of any failing elements are

smaller than the conditions given in Equations (55) and (56).

In summary, Equations (8) and (9) mark the transition from evolving microdamage to failure to macroscopic cracking. Prior to failure initiation, Equation (30) is

used to calculate the microdamage reduced ISV Sr, and the failure ISVs Sf

I , Sm

I ,

and Sm

II remain zero. Equations (21) and (22) are used to calculate the degraded

transverse and shear moduli. Subsequent to matrix failure initiation, microdamage

growth is precluded, and Sr remains at S∗

r , the value of Sr when Equation (8) was

satisfied. The degeneration of the transverse and shear moduli, resulting from matrix

transverse and shear cracking, is calculated using Equations (50) and (51). Finally

if Equation (9) is satisfied, the axial modulus is calculated using Equations (49) as

fiber breakage evolves in the element. Once the material moduli have been calculated

using the appropriate evolution equations, the stresses can be updated accordingly

using Equations (20).

3 Mesh Objectivity

The theory outlined in Section 2 eliminates the mesh dependency that arises from

ill-posedness when the elements enter the post-peak softening regime by introducing

a characteristic length into the formulation. The total SERR dissipated during the

evolution of the discontinuity is equal to the prescribed fracture toughness and is

independent of the element size. This approach, commonly referred to as the smeared

NASA/TM—2011-217401 21
crack approach, or crack band model, has been used to alleviate mesh dependency in

FEM since it was first developed by Baˇzant (1982) for post-peak strain softening in

concrete.

To exhibit the mesh objectivity of EST, a simple example is presented in this

section. One quarter of a 200 mm x 100 mm panel containing a central hole is

modeled with finite elements using the Abaqus, version 6.10-1 finite element software

Abaqus (2008). The panel contains a hole with a radius of rh = 5 mm in the center.

The left edge of the panel is constrained in the x-direction to simulate symmetry.

Similarly, the bottom edge is constrained in the y-direction. A uniform displacement

is applied to all the nodes on the top edge of the panel in y-direction. Details on

the panel geometry and boundary conditions are displayed in Figure 5. The panel is

composed of a generic, [90o], composite lamina with the fiber angle measured with

respect to the y-axis; thus, the applied displacement is perpendicular to the fiber

direction in the panel. EST is used to model damage and failure in the panel.

Four different meshes are used to evaluate the effect of mesh size on the response

of the panel. All four meshes consist of two-dimensional (2-D), plane stress, quadrilateral, S4R shell elements. The elements are linear, reduced integration elements and

contain four nodes and one integration point each. The density of the four meshes,

shown in Figure 6, increases within a region near the central hole. Average element

sizes equal 0.5rh, 0.2rh, 0.1rh, and 0.04rh are used in the four different meshes. Coarser

elements are used away from the hole to improve computation time. The four meshes

are subjected to the same boundary conditions and loading, and are composed of the

same material properties. The same elastic, damage and failure parameters are also

used in all four simulations.

The resultant, applied tensile stress (given by two times the sum of the reaction

forces at the nodes on the top edge divided by the cross-sectional area) normalized

by the transverse, mode I, critical strain times the transverse Young’s modulus  ̄σ is

NASA/TM—2011-217401 22
plotted versus the applied displacement normalized by the radius of the hole Δ for  ̄

the four different meshes in Figure 7. It can be seen that the mesh density has a

minimal effect on the load-deflection results. The small discrepancy in the results

between the four meshes can be attributed to the increased accuracy in the fields as

the mesh is refined. Moreover, the total energy dissipated is preserved from mesh to

mesh.

Figure 8 displays contours of the normalized, reduced microdamage ISV Sr immediately before failure initiation in the four different meshes. Sr is normalized by the

maximum Sr obtained in all four simulations which is 27% of the maximum allowable

Sr required to bring the moduli to zero. The four meshes exhibit similar microdamage

contours, but as the mesh is refined, the microdamage is contained to the vicinity of

the hole. Additionally, the global stress at which failure initiates reduces slightly as

the mesh is refined; this supports the previous statement that the discrepancies in

the stress-displacement responses were a facet of increasing field accuracy with mesh

refinement.

The transverse matrix failure degradation parameter Dm

I is plotted for all four

meshes in Figure 9. Figures 9a-9d show the failure pattern at the ultimate, global

load. The coarsest mesh shows that a crack band has grown at the intersection of

the hole and the bottom symmetric boundary and is propagating towards the free

edge, while moving away from the bottom boundary when the specimen ultimate

load is achieved. In the finer meshes, the crack band path is different. In Figures

9b-9c the crack band initiates at the hole slightly above the bottom boundary and

propagates towards the free edge while approaching the bottom boundary. The finest

mesh, Figure 9d, exhibits multiple crack bands near the hole, but only one of the

crack bands grows significantly. While in all simulations the first crack band initiates

at the hole, Figures 9c and 9d show some crack bands beginning to initiate at the free

edge and propagate inwards by the time the ultimate load is attained. Figures 9e-9h

NASA/TM—2011-217401 23
display the transverse crack band pattern once the specimen has lost all load-carrying

capability. The solution for the simulation with the finest mesh 0.04rh diverged before

all the load carrying capability was lost; so, Figure 9h presents the crack band path

at the final converged state, which is still far below the ultimate load state. In

Figures 9e-9h the same crack band patterns that developed in Figures 9a-9d are

evident, except those crack bands have saturated to maximum degradation: Dm

I = 1.

Additionally, the crack bands advancing from the free edge observed in Figures 9c and

9d have progressed further. However, this is well after the specimens reached their

ultimate load; therefore, it is assumed that those crack bands did not influence the

load carrying capability of the panels. The discrepancy in crack band path observed

for the different meshes indicate that the crack band path is dictated by the accuracy

of the fields surrounding the leading boundary of the crack band path.

4 Example - Center Notched Panels Subjected to

Uniaxial Tension

4.1 Experimental Details

Two center-notched panel (CNP) configurations were tested at the NASA Langley

Research Center (LaRC) [Bogert et al. (2006); Satyanarayana et al. (2007)]. The

geometrical details of the panel and testing boundary conditions are presented in

Figure 10. The panels were 3” wide and 11.5” long. Two 3” x 2.75” tabs were placed

on both ends of the specimens, leaving a gage section of 3” x 6” which is displayed

in Figure 10. A central notch was machined in each panel that was 0.75” wide and

had a notch tip radius of 0.09375”. The end tabs were clamped and a vertical,

tensile displacement (in the y-direction) was applied to the top tab using a servohydraulic testing machine. The bottom tab was fixed preventing any y-displacement

NASA/TM—2011-217401 24
of the bottom boundary of the gage section. The gripped tabs also prevented any

displacement in the x-direction at the top and bottom boundaries of the gage section.

The panels were comprised of laminated T800/3900-2 carbon fiber/toughened

epoxy composites. Three different lay-up configurations were tested; however one of

the configurations exhibited significant delamination. Since the focus of this work

is modeling in-plane damage and failure mechanisms, this configuration is not considered here. The two remaining configurations are presented in Table 1. The first

lay-up, Laminate-1, consists of 12 0◦ plies, and the second, Laminate-2, is a symmetric, multi-angle lay-up with 40% |45◦|, 40% 0◦, and 20% 90◦ layers.

Several strain gages where affixed to the test panel, labeled Sg-1 through Sg-4 in

Figure 10. Sg-1 was placed in the center of the panel, 1.5” above the notch. Sg-2 was

placed 1.5” above the notch tip. Sg-3 was attached in front of the notch, 0.5” from the

free edge, and Sg-4 placed at the notch tip. Global load versus displacement data, and

local strain gage data was reported in Refs. [Bogert et al. (2006); Satyanarayana et al.

(2007)], along with a post-test C-Scan of Laminate-1 and photograph of Laminate-1

and photograph of Laminate-2.

4.2 Finite Element Model Details

The linear elastic properties of T800/3900-2 used in the FEM models are presented

in Table 2, and were taken from Ref. [Bogert et al. (2006)]. The shear microdamage function gs utilized in Equation (22) was obtained from [45◦/-45◦]3S angle-ply

T800/3900-2 coupon tests. The transverse, tensile and compressive microdamage

functions were inferred by scaling the coefficients of the microdamage curves presented by Sicking (1992) for AS4/3502 by the ratio of the virgin transverse modulus

of T800/3900-2 to that of AS4/3502, as the stress-strain curves of the coupon laminates necessary to characterize es were not available. The polynomial forms of es and

NASA/TM—2011-217401 25
gs are

es(Sr) = es0 + es1Sr + es2S2

r + es3S3r + es4S4r + es5S5r (57)

gs(Sr) = gs0 + gs1Sr + gs2S2

r + gs3S3r + gs4S4r + gs5S5r (58)

The coefficients of the microdamage curves are presented in Table 3, and the curves

are plotted in Figure 12.

To increase computational efficiency, the first derivatives of the higher order microdamage polynomials are approximated using linear spline interpolants. Thus, the

microdamage evolution equation, Equation (30), is always second order in Sr yielding

a very efficient analytical solution. Since the value of Sr from the previous increment

is used to estimate which spline regime should be used to solve for Sr in the current

increment, the solution is checked to ensure that it falls within the applicable range

of the spline that was used. If it falls outside of the range of Sr that are valid for

the splines, the solution is calculated again using splines that accord to the solution

of the previous iteration. This procedure continues until the solution of Equation

(30) falls within the relevant range of Sr for the splines used in Equation (30). A

maximum number of iterations can be set, after which Equation (30) is solved using

the full polynomial forms of the damage functions, given in Equations (57) and (58),

by finding the eigenvalues of the companion matrix of the polynomial coefficients of

the evolution equation.

The axial mode I, transverse mode I, and shear mode II critical cohesive strains,

and fracture toughness’ are given in Table 4. The matrix mode I and mode II cohesive

critical strains (YT , YC, and Z) and the fracture toughnesses (Gm

IC and GmIIC) were

calibrated using data from Laminate-1. These values were adjusted until the global

load versus local strain at strain gage Sg-1 (see Figure 10) obtained from the model

provided the best match to the experimental data. Both the simulation results and

experimental data are presented in Figure 14a. Laminate-1 did not exhibit any axial

NASA/TM—2011-217401 26
failure; so, the fiber mode I parameters (XT and Gf

IC) were calibrated such that

the ultimate load from the simulation of Laminate-2 corresponded with the ultimate

load reported by Bogert et al. (2006) for Laminate-2. Subsequent work will outline a

procedure for measuring these values experimentally.

Displacement was applied to both laminates using the *DYNAMIC keyword in

Abaqus with the parameter APPLICATION = QUASI-STATIC. This implicit dynamic

solver is recommended for quasi-static problems exhibiting a high-degree of nonlinearity. This procedure uses numerical damping to stabilize the problem. The numerical

damping does not significantly affect the simulation results because the velocities in

these simulations are low. For Laminate-1 a total displacement of 0.0236”, and for

Laminate-2 a displacement of 0.0472”, is applied over 1000 seconds. The panels were

assigned a representative density of 0.057 lb/in.3. This implicit, dynamic technique

has advantages over traditional static, implicit solvers which have difficulty converging

when the material exhibits post-peak softening [Belytschko et al. (2000); Belytschko

and Mish (2001)], and is not limited by a minimum stable time step required with

explicit solvers [Hughes (2000)].

4.3 Results - Laminate-1

Global load P versus displacement Δ of a 4” section of Laminate-1 is compared

to results from the EST simulation in Figure 13. Very good agreement between the

model and the experimental results is achieved. The response of the specimen appears

to be linear until near 8,000 lbf., where the specimen begins deforming non-linearly.

The EST simulation captures the initiation and progression of the global nonlinearity

accurately. This panel was not loaded until catastrophic failure; hence, the data

presented in Figure 13 represents load versus displacement data prior to the ultimate

load of the specimen.

Local strain gage data (global load P versus local y-direction strain yy) from

NASA/TM—2011-217401 27
Laminate-1 is plotted with the results from the EST FEM model in Figure 14; please

refer to Figure 10 for locations of strain gages. Strain relaxation is observed in the gage

farthest away from the notch: Sg-1, shown in Figure 14a. The mode I and mode II

matrix failure parameters in EST were calibrated such that the model demonstrates

the same transition into strain relaxation at this location and at a similar global

applied load. This load, taken as the splitting load, is 8,250 lbf. in experiment and

8,210 lbf. in the model (summarized in Table 5). The transition to strain relaxation

is more abrupt in the experiment as evidenced by the sharp knee in the load-strain

curve, whereas, the transition in the model is more gradual. Prior to strain relaxation

at this point, the experiment displayed slight stiffening not observed in the model.

Additionally, the model response is much smoother than that of the experiment in

the strain relaxation regime. Even though the global loading is quasi-static, local

events, such as cracking, may be dynamic; therefore, the discrepancy in the strain

relaxation portion of the load-strain curves could be a result of dynamic matrix crack

growth and arrest in the test specimen. Local crack dynamics were not taken into

account in the model. Additionally, the jaggedness of the experimental data may be

a facet of the stochastics related to the local microstructure of the composite that

are not included in the model. The data from the experiment and simulation for

Sg-2, which is located 1.5” directly above the notch, are presented in Figure 14b.

The model predicts less strain at Sg-2, for a given load, than the experiment, but the

non-linear trends are very similar. This gage lies directly in front of the splitting crack

path, shown in Figure 10, and it is not realistic to expect perfect agreement in areas

experiencing high levels of damage and failure because of idealizations used to model

the evolution of cracks in the simulations. Figure 14c displays data for Sg-3, located in

front of the notch near the free edge. Very good agreement between the experimental

and simulation results are exhibited. The model accurately captures the non-linear

evolution of strain, away from the highly damaged regions, as a function of applied

NASA/TM—2011-217401 28
load. Finally, results for Sg-4 (located directly at the notch tip) are given in Figure

14d and includes the experiment and simulation display of axial strain relaxation.

As with Sg-2, Sg-4 shows less strain for a given applied load. However, the load at

which the strain at Sg-4 relaxes in both the experiment and model correlate well, in

accordance with the splitting load. Again, the relaxation response of the experiment

is discontinuous, but the model exhibits continuous behavior.

A C-Scan of the failed Laminate-1 specimen is displayed in Figure 15. Four

splitting cracks can be observed propagating outward from the notch tip, parallel to

the loading direction, towards the gripped edges. Contour plots of the normalized

microdamage obtained from the simulation are presented in Figure 16 at the splitting

load 8,200 lbf. and at 16,400 lbf. In these plots, Sr is normalized by the maximum

achievable value, Smax

r = 7.57 psi 1

3 , obtained from Figure 12. In Laminate-1, Sr

reached a maximum value equal to 0.171Smax

r . At the splitting load, the regions of

maximum damage are localized to small regions, along the same crack path observed

in Figure 15, embedded in a more widespread domain and exhibiting less severe

microdamage. A similar microdamage contour is observed at P = 16,400 lbf., except

the localized damage region has nearly proceeded to the fixed boundaries of the panel.

Figure 17 shows the shear failure degradation factor Dm

II at the splitting load and

16,400 lbf. The shear failure localizes into crack bands that are a single element wide

and progress equivalently to the cracks observed in the experiment. At the splitting

load (Figure 17a), the crack bands have progressed less than half of the way between

the notch and the panel boundary on either side of the notch. In Figure 17b, the

crack bands have nearly reached the fixed grip boundaries. Additionally, Figure 17b

displays some irregularity in the crack path. In these regions, the mesh is not uniform

because it needed to accommodate larger elements used to represent the strain gages

(see Figure 11). No axial failure was observed in this simulation. The shear crack path

was in a state of transverse compressions, and according to the traction-separation

NASA/TM—2011-217401 29
relations used in Figure 4b, transverse failure does not progress under compressive

conditions. Thus, contours of Df

I and Dm

I are not shown.

4.4 Results - Laminate-2

Numerical results for applied load versus displacement of a 4” section of Laminate-2

are presented in Figure 18. The experimental ultimate load 15,300 lbf. correlates well

(axial failure parameters were calibrated to obtain an ultimate load that most closely

matched the experimental data) with the ultimate load obtained from the model, also

15,300 lbf., and is summarized with the splitting load from the Laminate-1 analysis in

Table 5. The global response up to failure is nearly linear and failure occurs suddenly

and catastrophically.

Figure 19 compares the applied load versus strain gage results from the model

to the data from the experiment. Sg-1 and Sg-2 exhibited similar behavior; the

strain increases until the ultimate load is obtained, after which the strain relaxes

abruptly. The experimental data and numerical results both display this behavior.

The model exhibits slightly more strain, for a given load, prior to ultimate failure.

At Sg-3, the model predicts strain localization after the ultimate load is achieved.

The gage data shows a slight reduction in strain as the load drops; however, the gage

was placed directly in the crack path and may have been damaged when the panel

failed. The model results and experimental data for Sg-4 exhibit similar trends, but

the strain gage shows a large degree of nonlinearity at the notch tip. Bogert et al.

(2006) attributed this observed nonlinearity to local interlaminar stresses near the

notch free edge which caused some local delaminations. Since the focus of this work

was modeling in-plane damage mechanisms, these effects are not captured; however,

the model could be easily extended to incorporate delamination by placing DCZM

elements between continuum shell layers [Satyanarayana et al. (2007)].

In the experiment, the gages measure the strain over a continuous area associated

NASA/TM—2011-217401 30
with the size of the gage, but in the model, the strain is taken at the integration point

of an element; thus, these measures should not be expected to correspond exactly.

In areas where there are large gradients present, such as near a notch tip (Sg-4) or

near cracks (Laminate-1, Sg-2 or Laminate-2, Sg-3), it becomes even more difficult to

relate the strain gage data to numerical strains from a discretized continuum. This

may contribute to some of the discrepancies between the local gage data and the

model results in Figures 14 and 19.

A photograph taken of the failed, Laminate-2 specimen is presented in Figure 20.

The photograph shows that two macroscopic cracks initially propagate from the notch

tip towards the free edges, perpendicular to the applied load, in a self-similar fashion.

Eventually, the cracks turn and proceed towards the free edge at an angle. Bogert

et al. (2006) claim, supported by visual image correlation displacement data, that

there was some eccentricity in the specimen alignment, which resulted in deviation

from self-similar crack growth.

Normalized microdamage contours just prior to the ultimate load are presented for

the outermost 45◦, 0◦, -45◦, and 90◦ plies in Figure 21. Similar microdamage patterns

are evident in the 45◦ and -45◦ layers. Microdamage propagates outward, toward the

free edge, from the notch tip in petal-like patterns. The microdamage in these layers

is highly distributed throughout the plies. The 0◦ ply displays a more contained

microdamage pattern along the lines of the microdamage contours associated with

axial splitting as shown in Figure 16. A moderate level of microdamage is also

displayed in the 90◦ layer, but a low degree of microdamage is distributed throughout

most of the layer.

Figure 22 shows the axial failure degradation parameter Df

I at the ultimate load for

the four unique layers. A small amount of axial failure in the 45◦, 0◦, and -45◦ layers

can be observed at the notch tips. It appears that more failure occurs at one notch

tip than the other. This can be attributed to numerical imperfections resulting from

NASA/TM—2011-217401 31
dissimilar meshes at the opposite notch tips, that is the mesh is not symmetric about

the y-axis. No axial failure is observed in the 90◦ layer. Contours of the transverse,

mode I, failure degradation parameter Dm

I at the ultimate load are plotted in Figure

23. The failure patterns are similar in the 45◦ and -45◦ plies in Figures 23a and 23c

and are comparable to the microdamage contours in Figures 21a and 21c, except

the failure is restricted to regions on either side of the notch. Furthermore, small,

highly degraded domains can be observed propagating from the notch tip at an angle

corresponding to the fiber direction in the ply. The 90◦ layer exhibits some moderate

degradation in a localized region around the notch tips, and the 0◦ layer does not

exhibit much Dm

I . Contours of the shear, mode II, failure degradation parameter DmII

are presented at the ultimate load in Figure 24. Very similar failure paths can be

seen in the 45◦ and -45◦ layers and the patterns are nearly symmetric across both

centerlines of the panel. This is expected because as Figure 4c indicates, the sign of

the local shear strain does not affect the failure degradation. Dm

II in the 0◦ and 90◦

is limited to very small regions surrounding the notch tips.

Contours representing the microdamage in the four unique layers are presented

in Figure 25 after the panel has completely failed and lost all of its load carrying

capability. Although further matrix microdamage evolution is prohibited in elements

that have failed (transverse/shear or axial), in the other elements that have not failed,

matrix microdamage evolution continues. Nearly the entire 45◦ and -45◦ layers reach

a microdamage level of 0.18Smax

r . The 0◦ and 90◦ plies exhibit similar microdamage

patterns; however, low levels of microdamage are more widespread in the 90◦ ply.

Figure 26 shows the fiber failure path once the specimen has completely failed. All of

the layers, except the 90◦ layer, show self similar cracks propagating from the notch

tips towards the free edges of the panel. The angled crack path shown in Figure

20 was not reproduced because the eccentric loading (suspected in the test) was not

introduced into the simulation; therefore, the crack growth remained self-similar. A

NASA/TM—2011-217401 32
high degree of transverse matrix failure can be seen in the axial crack path in the

45◦, -45◦, and 90◦ plies in Figure 27. In the 0◦ layer, some transverse failure is

observed surrounding the fiber failure, as well as away from the axial failure path,

which resulted from a stress wave reflecting off the free edges when the axial crack

band reaches the boundary. Finally, Dm

II is presented after the specimen has failed

in Figure 28. Similar failure to Figures 24a and 24c in the 45◦ and -45◦ is exhibited,

but a highly degraded region has localized in the axial crack path. Figures 28b and

28d show fairly extensive regions containing a high degree of shear matrix failure

surrounding the axial failure path.

5 Conclusions

A thermodynamically-based, work potential theory for damage and failure in composite materials, enhanced Schapery theory (EST), was developed. A marked distinction

between damage and failure was introduced. Damage was considered to be the evolution of mechanisms that cause structural changes in the material such that the

non-linear tangent stiffness tensor remains positive definite. Failure was taken to be

the effect of structural changes in the material that result in loss of positive definiteness of the tangent stiffness matrix and post-peak strain softening. Separate internal

state variables (ISVs) were used to account for damage and three in-plane failure

mechanisms.

In EST, matrix microdamage, which includes matrix microcracking, shear banding, and microvoid growth, is responsible for all damage in a composite lamina and

was accounted for with a single ISV, along the lines of the original Schapery theory (ST) formulation. The relationship between the transverse and shear moduli

of the lamina were related to the ISV through a pair of experimentally-obtainable

microdamage functions.

NASA/TM—2011-217401 33
Three major, in-plane failure mechanisms applicable to continuous fiber-reinforced,

laminated, polymer matrix composites were identified: mode I matrix cracks, mode

II matrix cracks, and fiber breakage. A failure initiation criterion was used to mark

the transition from a damaging continuum to a damaged continuum with an embedded discontinuity. After failure initiation, microdamage evolution ceases and separate ISVs are introduced to incorporate the effects of the three major failure mechanisms. Evolution of the failure ISVs is based upon traction-separation laws (which

are a functions of the appropriate fracture toughnesses) and a characteristic element

length. Typically, the existence of a non-positive definite stiffness tensor would result in pathologically mesh dependent solutions; however, in EST, mesh objectivity

is ensured by incorporating a characteristic length scale into the failure evolution.

In Section 3 of this paper, mesh objectivity is demonstrated. A unidirectional

composite plate with a central hole obeying EST is loaded in transverse tension and

the response is calculated using four different, increasingly refined, meshes. The global

stress versus strain response remained unaffected by the change in mesh, save for the

effects from increased accuracy of local fields in the vicinity of the hole with denser

meshing.

Two center-notched panels, with different lay-ups, composed of T800/3900-2 were

tested under tensile loading at NASA LaRC. Global load versus displacement and

global load versus local strain gage strain data were compared to results obtained

from FEM models utilizing EST in Section 4. Quantitatively, very good correlation

was achieved for both laminates. Furthermore, damage and failure paths predicted

by the models matched well with the experimental results.

NASA/TM—2011-217401 34
ID Stacking Sequence Thickness (in.)

Laminate-1 [0◦]12 0.078

Laminate-2 [45◦/0◦/-45◦/0◦/90◦]S 0.065

Table 1: T800/3900-2 lay-up configurations used in CNP tests at NASA LaRC.

Property T800/3900-2

E11 (Msi) 23.2

E22 (Msi) 1.3

G12 (Msi) 0.9

ν12 0.28

Table 2: Linear elastic properties for T800/3900-2 used in FEM models.

es(Sr) Coefficients Tension Compression gs(Sr) Coefficients

es0 1.0000 1.0000 gs0 1.0000

es1 -6.0373E-2 8.4887E-4 gs1 -3.0567E-2

es2 2.5937E-2 2.8002E-2 gs2 -1.2135E-1

es3 -1.5789E-2 -6.2122E-3 gs3 3.7438E-2

es4 2.2571E-3 N/A gs4 -4.5405E-4

es5 -1.0440E-4 N/A gs5 1.9532E-4

Table 3: Microdamage function coefficients for T800/3900-2 used in FEM models.

Property T800/3900-2

XT 0.021

YT 0.0048

YC 0.0119

Z 0.0075

Gf

IC 1026 lbf.-in.

in.2

Gm

IC 170.0 lbf.-in.

in.2

Gm

IIC 13.54 lbf.-in.

in.2

Table 4: Failure parameters for T800/3900-2.

Type Experimental Numerical

Laminate-1 Splitting 8,250 lbf. 8,210 lbf.

Laminate-2 Ultimate 15,300 lbf. 15,300 lbf.

Table 5: Critical experimental and simulation loads for Laminate-1 and Laminate-2.

NASA/TM—2011-217401 35
Figure 1: Typical stress-strain curve, containing pre-peak nonlinearity and post-peak

strain softening, showing the total elastic (WE) and total dissipated (WS) potentials.

Figure 2: Typical stress-strain curve with a positive-definite tangent stiffness exhibiting microdamage, showing the elastic (W) and irrecoverable (S) portions.

NASA/TM—2011-217401 36
Figure 3: Schematic showing the transition form a continuum to a cohesive zone due

to the initiation of macroscopic cracks. The essential, constitutive variables switch

from stress and strain to traction and separation.

(a) Mode I fiber fracture. (b) Mode I matrix fracture.

(c) Mode II matrix fracture.

Figure 4: Triangular traction versus separation which dictates the behavior of cohesive cracks embedded in the continuum. The total area under the traction-separation

law represents the material fracture toughness Gj

mC. The area above the unloading

line for a given traction-separation state is the strain energy release rate Gj

m.

NASA/TM—2011-217401 37
Figure 5: Example problem used to demonstrate mesh objectivity of EST. One

quarter of a 200 mm x 100 mm containing a central hole with a radius of 5 mm

is loaded in tension with symmetric boundary conditions on the bottom and left

boundaries.

(a) 0.5rh. (b) 0.2rh. (c) 0.1rh. (d) 0.04rh.

Figure 6: Four mesh densities used to demonstrate the mesh objectivity of EST.

NASA/TM—2011-217401 38
0 0.05 0.1 0.15 0.2 0.25 0.3 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Δ ̄

 ̄σ

0.5rh

0.2rh

0.1rh

0.04rh

Figure 7: Reaction stress normalized by critical axial strain times axial Young’s

modulus verse applied displacement normalized by hole radius for four different mesh

densities.

(a) 0.5rh,  ̄σ = 0.50. (b) 0.2rh,  ̄σ =

0.42.

(c) 0.1rh,  ̄σ =

0.37.

(d) 0.04rh,  ̄σ =

0.36.

Figure 8: Contours of the reduced microdamage ISV Sr, normalized by the maximum Sr obtained from all simulations, immediately prior to failure initiation for four

different mesh densities.

NASA/TM—2011-217401 39
(a) 0.5rh,  ̄σ = 0.80. (b) 0.2rh,  ̄σ = 0.81.

(c) 0.1rh,  ̄σ = 0.82. (d) 0.04rh,  ̄σ = 0.82.

(e) 0.5rh,  ̄σ = 0.030. (f) 0.2rh,  ̄σ = 0.018.

(g) 0.1rh,  ̄σ = 0.042. (h) 0.04rh,  ̄σ = 0.242.

Figure 9: Contours of the transverse degradation parameter Dm

I , indicative of the

transverse crack path in the specimens. The contours (a-d) are presented at the

ultimate load achieved by the specimens and (e-h) after the specimens have lost their

load carrying capability.

NASA/TM—2011-217401 40
Figure 10: Geometry, boundary conditions, and strain gage (Sg) locations of CNPs

tested at NASA LaRC [Bogert et al. (2006)].

Figure 11: FEM mesh used to simulate tensile loading of CNPs.

NASA/TM—2011-217401 41
0 1 2 3 4 5 6 7 8 0

0.2

0.4

0.6

0.8

1

1.2

Sr (psi 1

3 )

gs

(a) Shear microdamage function obtained

from [±45◦]S angle-ply laminate.

0 1 2 3 4 5 6 7 8 0

0.2

0.4

0.6

0.8

1

1.2

Sr (psi 1

3 )

es

Tension

Compression

(b) Transverse tension and compression

microdamage functions obtained by scaling data for AS4/3502 in Ref. [Sicking

(1992)].

Figure 12: Microdamage functions for T800/3900-2 used in FEM models.

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014

0

2000

4000

6000

8000

10000

12000

14000

Δ of a 4” Section (in.)

P (lbf.)

Experiment

EST

Figure 13: Applied load versus displacement of a 4” section for Laminate-1.

NASA/TM—2011-217401 42
0 1 2 3 4 5 6 7 8

x 10−4

0

2000

4000

6000

8000

10000

12000

14000

yy

P (lbf.)

Experiment

EST

(a) Sg-1.

0 0.5 1 1.5 2 2.5

x 10−3

0

2000

4000

6000

8000

10000

12000

14000

yy

P (lbf.)

Experiment

EST

(b) Sg-2.

0 0.5 1 1.5 2 2.5

x 10−3

0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

yy

P (lbf.)

Experiment

EST

(c) Sg-3.

0 1 2 3 4 5

x 10−3

0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

yy

P (lbf.)

Experiment

EST

(d) Sg-4.

Figure 14: Applied load versus local strain for Laminate-1.

Figure 15: C-Scan of failed Laminate-1 specimen [Bogert et al. (2006)].

NASA/TM—2011-217401 43
(a) P = 8,210 lbf.

(b) P = 16,400 lbf.

Figure 16: Normalized matrix microdamage contour Sr

Smax

r

in Laminate-1.

NASA/TM—2011-217401 44
(a) P = 8,210 lbf.

(b) P = 16,400 lbf.

Figure 17: Matrix shear failure degradation Dm

II in Laminate-1.

NASA/TM—2011-217401 45
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035

0

2000

4000

6000

8000

10000

12000

14000

16000

Δ of a 4” Section (in.)

P (lbf.)

Experiment

EST

Figure 18: Applied load versus displacement of a 4” section for Laminate-2.

0 1 2 3 4 5 6

x 10−3

0

2000

4000

6000

8000

10000

12000

14000

16000

yy

P (lbf.)

Experiment

EST

(a) Sg-1.

0 1 2 3 4 5 6 7

x 10−3

0

2000

4000

6000

8000

10000

12000

14000

16000

yy

P (lbf.)

Experiment

EST

(b) Sg-2.

0 0.002 0.004 0.006 0.008 0.01 0.012 0

2000

4000

6000

8000

10000

12000

14000

16000

yy

P (lbf.)

Experiment

EST

(c) Sg-3.

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0

2000

4000

6000

8000

10000

12000

14000

16000

yy

P (lbf.)

Experiment

EST

(d) Sg-4.

Figure 19: Applied load versus local strain for Laminate-2.

NASA/TM—2011-217401 46
Figure 20: Photograph of failed Laminate-2 specimen [Bogert et al. (2006)].

(a) 45◦ Layer. (b) 0◦ Layer.

(c) -45◦ Layer. (d) 90◦ Layer.

Figure 21: Normalized matrix microdamage contour Sr

Smax

r

in Laminate-2 just prior

to first axial failure initiation P = 8,640 lbf.

NASA/TM—2011-217401 47
(a) 45◦ Layer. (b) 0◦ Layer.

(c) -45◦ Layer. (d) 90◦ Layer.

Figure 22: Fiber failure degradation Df

I in Laminate-2 at ultimate load P = 15,300

lbf (magnified view of region near notch).

NASA/TM—2011-217401 48
(a) 45◦ Layer. (b) 0◦ Layer.

(c) -45◦ Layer. (d) 90◦ Layer.

Figure 23: Transverse matrix failure degradation Dm

I in Laminate-2 at ultimate load

P = 15,300 lbf.

NASA/TM—2011-217401 49
(a) 45◦ Layer. (b) 0◦ Layer.

(c) -45◦ Layer. (d) 90◦ Layer.

Figure 24: Shear matrix failure degradation Dm

II in Laminate-2 at ultimate load P

= 15,300 lbf.

NASA/TM—2011-217401 50
(a) 45◦ Layer. (b) 0◦ Layer.

(c) -45◦ Layer. (d) 90◦ Layer.

Figure 25: Normalized matrix microdamage contour Sr

Smax

r

in Laminate-2 after specimen has lost load carrying capability.

NASA/TM—2011-217401 51
(a) 45◦ Layer. (b) 0◦ Layer.

(c) -45◦ Layer. (d) 90◦ Layer.

Figure 26: Fiber failure degradation Df

I in Laminate-2 after specimen has lost load

carrying capability.

NASA/TM—2011-217401 52
(a) 45◦ Layer. (b) 0◦ Layer.

(c) -45◦ Layer. (d) 90◦ Layer.

Figure 27: Transverse matrix failure degradation Dm

I in Laminate-2 after specimen

has lost load carrying capability.

NASA/TM—2011-217401 53
(a) 45◦ Layer. (b) 0◦ Layer.

(c) -45◦ Layer. (d) 90◦ Layer.

Figure 28: Shear matrix failure degradation Dm

II in Laminate-2 after specimen has

lost load carrying capability.

NASA/TM—2011-217401 54
References

Abaqus (2008), Abaqus User’s Manual, Vol. 1-3, Version 6.10-1, Dassault Syst`emes

Simulia Corp., Providence, RI.

Allen, D. H., C. E. Harris, and S. E. Groves (1987a), A thermomechanical constitutive

theory for elastic composites with distributed damage - i. theoretical development,

Int. J. Solids Struct., 23 (9), 1301–1318.

Allen, D. H., C. E. Harris, and S. E. Groves (1987b), A thermomechanical constitutive

theory for elastic composites with distributed damage - ii. aplication to matrix

cracking in laminated composites, Int. J. Solids Struct., 23 (9), 1319–1338.

Basu, S. (2005), Computational modeling of progrssive failure and damage in composite laminates, Ph.D. thesis, University of Michigan, Ann Arbor, MI.

Basu, S., A. M. Waas, and D. R. Ambur (2006), Compressive failure of fiber composites under multiaxial loading, Int. J. Solids Structures, 44 (9), 2648–2676.

Baˇzant, Z., and L. Cedolin (1979), Blunt crack band propagation in finite element

analysis, J. Eng. Mech. Div.-ASCE, 105, 297–315.

Baˇzant, Z. P. (1982), Crack band model for fracture of geomaterials, in Proceedings of the 4th International Conference on Numerical Methods in Geomechanics,

Edmonton, Canada.

Baˇzant, Z. P., and B. H. Oh (1983), Crack band theory for fracture of concrete, Mater.

and Struct., 16, 155–77.

Bednarcyk, B. A., J. Aboudi, and S. M. Arnold (2010), Micromechanics modeling of

composites subjected to multiaxial progressive damage in the constituents, AIAA

J., 48, 1367–1378.

Belytschko, T., and K. Mish (2001), Computability in non-linear solid mechanics, Int.

J. Numer. Methods, 52, 3–21.

Belytschko, T., W. K. Liu, and B. Moran (2000), Nonlinear finite elements for continua and structures, Wiley, New York.

Bogert, P. B., A. Satyanarayana, and P. B. Chunchu (2006), Comparison of damage

path predicions for composite laminates by explicit and standard finite element

analysis tool, in 47th AIAA Structures, Structural Dynamics, and Materials Conference.

de Borst, R., and P. Nauta (1985), Non-orthogonal cracks in a smeared finite element

model, Eng. Comput., 2, 35–46.

Dvorak, G. J., N. Laws, and M. Hejazi (1985), Analysis of progressive matrix cracking

in composite laminates i. thermoelastic properties of a ply with cracks, J. Compos.

Mater., 19, 216–234.

NASA/TM—2011-217401 55
Green, B. G., M. R. Wisnom, and S. R. Hallett (2007), An experimental investigation

into the tensile strength scaling of notched composites, Compos.: Part A, 38, 867–

878.

Gudmundson, P., and S. Ostlund (1992), First order analysis of stiffness reduction  ̈

due to matrix cracking, J. Compos. Mater., 26, 1009–1030.

Gustafson, P. A. (2008), Analytical and experimental methods for adhesively bonded

joints subjected to high temperatures, Ph.D. thesis, University of Michigan, Ann

Arbor, MI.

Gustafson, P. A., and A. M. Waas (2009), The influence of adhesive constitutive

parameters in cohesive zone finite element models of adhesively bonded joints, Int.

J. Solids Struct., pp. 2201–2215.

Hallett, S. R., W.-G. Jiang, B. Khan, and M. R. Winsom (2008), Modelling the interaction between matrix cracks and delamination damage in scaled quasi-isotropic

specimens, Compos. Sci. and Technol., 68, 80–89.

Hashin, Z., and A. Rotem (1973), A fatigue failure criterion for fiber reinforced composite materials, J. Composite Materials, 7, 448–464.

Hinterhoelzl, A., and R. A. Schapery (2004), Fem implementation of a threedimensional visoelastic constitutive model for particulate composites with damage

growth, Mech. Time-Depend. Mat., 8 (1), 65–94.

Hoek, E., and Z. T. Bieniawski (1965), Brittle rock fracture propagation in rock under

compression, Int. J. Fract., 1 (3), 137–155.

Hughes, T. J. R. (2000), The Finite Element Method: Linear Static and Dynamic

Finite Element Analysis, Dover Publications, Inc.

Kachanov, L. M. (1958), On the creep fracture time, Izvestiya Akademii Nauk SSSR,

Otdeleniya Tekhnika Nauk, 8 (1), 26–31.

Kachanov, L. M. (1986), Introduction to Continuum Damage Mechanics, MartinusNijhoff, Dordrecht, The Netherlands.

Krajcinovic, D. (1996), Damage Mechanics, Elsevier, New York.

Lamborn, M. J., and R. A. Schapery (1988), An investigation of deformation pathindependence of mechanical work in fiber-reinforced plastics, in Proceedings of the

Fourth Japan-U.S. Conference on Composite Materials, Technomic Publishing Co.,

Inc., Lancaster, PA.

Lamborn, M. J., and R. A. Schapery (1993), An investigation of the existence of a

work potential for fiber-reinforced plastic, J. Compos. Mater., 27, 352–382.

Laws, N., and G. J. Dvorak (1988), Progressive transverse cracking in composite

laminates, J. Compos. Mater., 22, 900–916.

NASA/TM—2011-217401 56
Lee, J.-W., D. H. Allen, and C. E. Harris (1989), Internal state variable approach for

predicting stiffness reductions in fibrous laminated composites with matrix cracks,

J. Compos. Mater., 23, 1273–1291.

Lemaitre, J. (1996), A Course on Damage Mechanics, 2nd Ed., Springer-Verlag,

Berlin, Heidelberg, New York.

Lemaitre, J., and J.-L. Chaboche (1994), Mechanics of Solid Materials, Cambridge

University Press.

Matzenmiller, A., J. Lubliner, and R. L. Taylor (1995), A constitutive model for

anisotropic damage in fiber-composites, Mech. Mater., 20 (2), 125–152.

McCartney, L. N. (1992a), Theory of stress transfer in a 0o-90o-0o cross-ply laminate

containing a parallel array of transverse cracks, J. Mech. Phys. Solids, 40 (1), 27–68.

McCartney, L. N. (1992b), Mechanics for the growth of bridged cracks in composite

materials: Part i. basic principles, J. Compos. Technol. & Res., 14 (3), 133–154.

McCartney, L. N. (1998), Predicting transverse crack formation in cross-ply laminates,

Compos. Sci. Technol., 58, 1069–1081.

Nairn, J. A. (1989), Strain energy release rate of composite microcracking: A variational approach, J. Compos. Mater., 23, 1106–1129.

Ng, W. H., A. G. Salvi, and A. M. Waas (2010), Characterization of the in-situ

non-linear shear response of laminated fiber-reinforced composites, Compos. Sci.

Technol., 70 (7), 1126–1134.

Noda, J., T. Okabe, N. Takeda, and M. Shimizu (2006), Tensile strength of cfrp crossply laminates containing transverse cracks, Adv. Comps. Mater., 15 (1), 81–93.

Paas, M. H. J. W., P. J. G. Schreurs, and W. A. M. Brekelmans (1992), A continuum

approach to brittle and fatigue damage: Theory and numerical procedures, Int. J.

Solids Struct., 30 (4), 579–599.

Pietruszczak, S., and Z. Mroz (1981), Finite element analysis of deformation of strainsoftening materials, Int. J. Numer. Methods Eng., 17, 327–334.

Pineda, E. J., A. M. Waas, B. A. Bednarcyk, C. S. Collier, and P. W. Yarrington

(2009), Progressive damage and faiure modeling in notched laminated fiber reinforced composites, Int. J. Fract., 158, 125–143.

Pineda, E. J., A. M. Waas, B. A. Bednarcyk, and C. S. Collier (2010), Computational

implementation of a thermodynamically based work potential model for progressive microdamage and transverse cracking in fiber-reinforced laminates, in 51st

AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials

Conference, Orlando, FL.

NASA/TM—2011-217401 57
Pinho, S. T., C. G. D ́avila, P. P. Camanho, and L. Iannucci (2005), Failure models

and criteria for frp under in-plane or three-dimensional stress states including shear

non-linearity, NASA/TM 2005-213530.

Puck, A., and H. Sch ̈urmann (1998), Failure analysis of frp laminates by means of

physically based phenomenological models, Comps. Sci. Technol., 58, 1045–1067.

Puck, A., and H. Sch ̈urmann (2002), Failure analysis of frp laminates by means of

physically based phenomenological models, Comps. Sci. Technol., 62, 1633–1622.

Rice, J. R. (1971), Inelastic constitutive relations for solids: an internal-variable

theory and its application to metal plasticity, J. Mech. Phys. Solids, 19, 433–455.

Satyanarayana, A., P. B. Bogert, and P. B. Chunchu (2007), The effect of delamination on damage path and failure load prediction for notched composite laminates,

in 48th AIAA Structures, Structural Dynamics, and Materials Conference.

Schapery, R. A. (1989), Mechanical characterization and analysis of inelastic composite laminates with growing damage, Mechanics & Materials Center Report 5762-

89-10, Texas A & M University, College Station, TX 77804.

Schapery, R. A. (1990), A theory of mechanical behaviour of elastic media with

growing damage and other changes in structure, J. Mech. Phys. Solids, 38 (2),

1725–1797.

Schapery, R. A. (1995), Prediction of compressive strength and kink bands in composites using a work potential, Int. J. Solids Structures, 32 (6), 739–765.

Schapery, R. A., and D. L. Sicking (1995), A theory of mechanical behaviour of

elastic media with growing damage and other changes in structure, in Mechanical

Behaviour of Materials, edited by A. Bakker, pp. 45–76, Delft University Press,

Delft, The Netherlands.

Sicking, D. L. (1992), Mechanical characterization of nonlinear laminated composites

with transverse crack growth, Ph.D. thesis, Texas A&M University, College Station,

TX.

Spencer, B. W. (2002), Finite elements with embedded discontinuities for modeling

reinforced concrete members, Ph.D. thesis, Brigham Young University, Provo, UT.

Talreja, R. (1985a), A continuum mechanics cahracterization of damage in composite

materials, P. Roy. Soc. Lond. A Mat., 4, 335–375.

Talreja, R. (1985b), Transverse cracking and stiffness reduction in composite laminates, J. Composite Materials, 19, 355–275.

Talreja, R. (Ed.) (1994), Composite Materials Series, vol. 9, Elsevier Science B.V.,

Amsterdam, The Netherlands.

NASA/TM—2011-217401 58
Tan, S. C., and R. J. Nuismer (1989), A theory for progressive matrix cracking in

composite laminates, J. Compos. Mater., 23, 1029–1047.

Turon, A., P. P. Camanho, J. Costa, and C. G. D ́avila (2006), A damage model for the

simulation of delamination in advanced composites under variable-mode loading,

Mech. Mater., 38 (11).

Voyiadjis, G., and P. I. Kattan (2005), Damage Mechanics, Taylor and Francis, Boca

Raton, FL.

Xie, D. E., A. Salvi, C. E. Sun, A. M. Waas, and A. Caliskan (2006), Discrete

cohesive zone model to simulate static fracture in 2-d triaxially braided carbon

fiber composites, J. Compos. Mater., 40 (22).

Yang, Q., and B. Cox (2005), Cohesive models for damage evolution in laminated

composites, Int. J. Frac., 133, 107–137.

NASA/TM—2011-217401 59
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

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-12-2011 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Numerical Implementation of a Multiple-ISV Thermodynamically-Based Work Potential 

Theory for Modeling Progressive Damage and Failure in Fiber-Reinforced Laminates 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Pineda, Evan, J.; Waas, Anthony, M. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 001098.04.02.03 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-18048 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2011-217401 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 24, 39, and 61 

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

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

A thermodynamically-based work potential theory for modeling progressive damage and failure in fiber-reinforced laminates is presented. The current, multiple-internal state 

variable (ISV) formulation, enhanced Schapery theory (EST), utilizes separate ISVs for modeling the effects of damage and failure. Damage is considered to be the effect of 

any structural changes in a material that manifest as pre-peak non-linearity in the stress versus strain response. Conversely, failure is taken to be the effect of the evolution of 

any mechanisms that results in post-peak strain softening. It is assumed that matrix microdamage is the dominant damage mechanism in continuous fiber-reinforced polymer 

matrix laminates, and its evolution is controlled with a single ISV. Three additional ISVs are introduced to account for failure due to mode I transverse cracking, mode II 

transverse cracking, and mode I axial failure. Typically, failure evolution (i.e., post-peak strain softening) results in pathologically mesh dependent solutions within a finite 

element method (FEM) setting. Therefore, consistent character element lengths are introduced into the formulation of the evolution of the three failure ISVs. Using the 

stationarity of the total work potential with respect to each ISV, a set of thermodynamically consistent evolution equations for the ISVs is derived. The theory is implemented 

into commercial FEM software. Objectivity of total energy dissipated during the failure process, with regards to refinements in the FEM mesh, is demonstrated. The model is 

also verified against experimental results from two laminated, T800/3900-2 panels containing a central notch and different fiber-orientation stacking sequences. Global load 

versus displacement, global load versus local strain gage data, and macroscopic failure paths obtained from the models are compared to the experiments.

15. SUBJECT TERMS

Fiber composites; Polymer matrix composites; Carbon fiber reinforced plastics; Computational mechanics; Damage; Fracture 

mechanics 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

66 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18