Evan J. Pineda and Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Craig S. Collier

Collier Research Corporation, Newport News, Virginia

Computational Implementation of a 

Thermodynamically Based Work Potential Model 

For Progressive Microdamage and Transverse 

Cracking in Fiber-Reinforced Laminates

NASA/TM—2012-217243

May 2012

AIAA–2010–2815
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 and Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Craig S. Collier

Collier Research Corporation, Newport News, Virginia

Computational Implementation of a 

Thermodynamically Based Work Potential Model 

For Progressive Microdamage and Transverse 

Cracking in Fiber-Reinforced Laminates

NASA/TM—2012-217243

May 2012

AIAA–2010–2815

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

51th Structures, Structural Dynamics, and Materials Conference

cosponsored by the AIAA, ASME, ASCE, AHS, and ASC

Orlando, Florida, April 12–15, 2010
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. 
Computational Implementation of a Thermodynamically Based

Work Potential Model For Progressive Microdamage and

Transverse Cracking in Fiber-Reinforced Laminates

A continuum-level, dual internal state variable, thermodynamically based, work potential model, Schapery

Theory, is used capture the effects of two matrix damage mechanisms in a fiber-reinforced laminated composite: microdamage and transverse cracking. Matrix microdamage accrues primarily in the form of shear

microcracks between the fibers of the composite. Whereas, larger transverse matrix cracks typically span the

thickness of a lamina and run parallel to the fibers. Schapery Theory uses the energy potential required to

advance structural changes, associated with the damage mechanisms, to govern damage growth through a set

of internal state variables. These state variables are used to quantify the stiffness degradation resulting from

damage growth. The transverse and shear stiffness’ of the lamina are related to the internal state variables

through a set of measurable damage functions. Additionally, the damage variables for a given strain state can

be calculated from a set of evolution equations. These evolution equations and damage functions are implemented into the finite element method and used to govern the constitutive response of the material points in

the model. Additionally, an axial failure criterion is included in the model. The response of a center-notched,

buffer strip-stiffened panel subjected to uniaxial tension is investigated and results are compared to experiment.

I. Introduction

D

AMAGE and failure analysis of composite materials continues to be prevalent and challenging research topic. The

intricacies of damage and failure in composite materials are attributed to numerous sources. The properties of

constituent materials, the architecture of the reinforcement material, the interfacial behavior, and the applied loading

configuration all influence the type and evolution of damage mechanisms in composite materials. Furthermore, many

globally observed damage mechanisms are a result of the interaction of the damage mechanisms in the different phases

and across various length scales. It is the interaction of damage mechanisms, due to the heterogeneity of the material,

that results in composites behaving more like a structure and is what separates them from monolithic materials. In

order to capture the non-linear behavior of composites correctly, the fundamental mechanisms must be modeled using

the underlying physics behind those mechanisms. Approaching composite modeling in this manner can lead to very

accurate, predictive models; which in turn, can lead to great weight and cost savings in design of composite structures.

For unidirectional, fiber-reinforced, epoxy laminates, numerous damage mechanisms are observed.1 Damage first

begins in the matrix of the composite at a scale on the order of the fibers. At low stress levels, matrix microdamage

arises in the form of micro-void growth, shear microcracks and shear bands between fibers.2, 3 This damage mechanism

is highly distributed and is the primary cause for the globally observed non-linear behavior up to the onset of more

severe damage mechanisms. As the load increases, transverse cracks begin to nucleate. These cracks typically run

parallel to the fibers and span the thickness of a lamina in the fiber-reinforced laminate (FRL).4–10 This mechanism

occurs on a larger length scale than microdamage. Additionally, the effects of this mechanism on the overall stiffness

1 of 22

Evan J. Pineda and Anthony M. Waas 

University of Michigan 

Ann Arbor, Michigan 48109 

Brett A. Bednarcyk 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Craig S. Collier 

Collier Research Corporation 

Newport News, Virginia 23607 

NASA/TM—2012-217243 1
of the composite are more pronounced. In angle-ply and unidirectional coupons, the initiation of transverse cracking

quickly leads to catastrophic failure of the specimen. Although these damage mechanisms are exclusive to the matrix

of the FRL, they are heavily influenced by the fiber properties and architectures. Placement of the fibers creates local

stress concentrations that nucleate and drive cracks within the materials. Also, the stiff fibers are responsible for

arresting matrix crack growth. Thus, the behavior of neat resin is different than the in-situ matrix

Axial (fiber-direction) damage of a ply in a laminate is attributed to damage in the fibers; however, the matrix

and fiber/matrix interface also play a role in axial damage. Carbon fibers, used in common FRL applications, usually

behave linearly until the fibers rupture suddenly. However in tension, axial failure of composite plies within a laminate

is often a combination of fiber breakage and the fiber being pulled out of the surrounding matrix.10–12 Fiber-pullout

can occur as macroscopic cracks progress within a laminate. It can result when there is a through-thickness matrix

crack, but only some, or none, of the fibers in that region have broken. Fiber-pullout can also arise when fiber

failure develops away from the primary fracture plane. In these scenarios the fibers must debond and “pull out” of

the surrounding matrix in order for the crack to open further. Thus, fiber-pullout is heavily dependent on the matrix

properties and interfacial properties. The onset of axial damage typically indicates that catastrophic failure of the

laminate is forthcoming; however the process can be slowed due to fiber-bridging between the two crack faces.

In compression, the globally observed axial damage mechanism is fiber micro-buckling.13–17 Small misalignments

of the fibers lead to fiber rotation as the structure is loaded. These rotations induce matrix softening in shear; which,

in turn, leads to further fiber rotations. Eventually, this leads to a runaway instability and a fiber kink band forms.

Often, in-plane lamina damage is accompanied by delamination between the layers. Interlamina damage occurs

when adjacent layers displace relative to each other and the bond between those layers is weakened. Delamination

can result from edge effects, or it may develop when a matrix crack, within a ply, reaches the interface of the adjacent

lamina. It has been shown experimentally that delamination can influence the globally observed failure modes of

open-hole tension, composite specimens.18–20 At the same time, in compression of open hole specimens, it has been

shown that fiber kink banding is dominant and leads to delamination formation.21–23

There exist numerous approaches to modeling damage mechanisms in a composite. Some methods for capturing

in-plane damage mechanisms at the continuum-level include using empirical failure criteria, continuum damage mechanics, and/or fracture mechanics and discontinuous methods.10, 12, 24–33 In-plane damage mechanisms can also be

captured with micromechanics models including repeating unit cell (RUC) and representative volume element (RVE)

techniques, as well as other methods.5, 34–39 Moreover, micromechanics models can be tied to the continuum scale

using multiscale methods in which the RUC is homogenized to give the effective response of the continuum.40–45

Commonly used techniques for introducing interlamina damage into a composite simulation include the virtual crack

closure technique (VCCT), cohesive zone model (CZM)10, 46–48 and discrete cohesive zone model (DCZM).49–52 Additionally, multiscale methods can also be used to simulate interlamina effects.53

The method presented in this work uses three different models to encapsulate the effects of various damage mechanisms in an FRL. Matrix microdamage and transverse cracking are modeled using a dual internal state variable (ISV)

formulation of the thermodynamically based work potential model developed by Schapery, and presented in Ref. 24.

Single and dual ISV formulations of Schapery Theory (ST) have been used to model microdamage alone and both

microdamage and transverse cracking in analytical and numerical frameworks.2, 17, 54, 55 Lamina-level ST represents

each layer as a homogenized continuum and uses the amount of energy necessary to advance structural changes within

the material, associated with the damage mechanisms, to govern damage growth in a thermodynamically consistent

manner. Each damage mechanism is accounted for through separate ISVs which are related to the degraded stiffness’

through a set of measurable damage functions. In this paper, ST is implemented into the finite element method (FEM)

and is used to control the constitutive behavior of the elements. Additionally, a simple maximum strain axial failure

criterion in used for axial failure. Upon satisfaction of the criterion, the stiffness of the corresponding material point

is significantly degraded. A more physically-based model for axial failure is desired and will be included in future

models. Finally, delamination is included in the model via cohesive elements. The details of the finite element model,

including the ST implementation, and axial failure are presented in Section IV.

Center-notched, IM7/8552 panels were tested in uniaxial tension at the NASA Langley Research Center.56 Two,

one inch sections, spanning the length of the panels, were reinforced with additional 0◦ plies. These sections, referred

to as buffer strips, were included to increase the ultimate load of the panels, without introducing too much additional

weight, and to change the failure mode of the panels. The experimental program is outlined in Section III. Section IV

provides details of an FEM model of the buffer strip panel. The global load versus displacement and local strain gage

data from the experiment are compared to the results of the FEM model in Section V. Additionally, analysis of the

damage and failure mechanisms observed in the FEM model is provided in Section V.

2 of 22

NASA/TM—2012-217243 2
II. Thermodynamically Based Work Potential Model

A thermodynamically based framework for modeling elastic solids undergoing structural changes was developed

by Schapery.24 In this theory the total work potential WT is partitioned into a recoverable, elastic portion W, and a

dissipated, or inelastic, part WS.

WT = W + WS (1)

As the material is loaded, a portion of the applied work potential facilitates structural changes in the material. The evolution of these structural changes affects the elastic properties of the material. However, a fraction of the total applied

work potential is recovered when the material is unloaded. The magnitude of energy recovered is contingent upon the

current, degraded elastic properties. It is assumed upon subsequent reloading, the material behaves linearly, exhibiting

the elastic properties observed during unloading, until it reaches the previous maximum strain state. After this state

is achieved, structural changes resume further degrading the elastic moduli of the material. This secant assumption

is not mandatory and may be relaxed depending on the material, but it is acceptable for fiber reinforced composites

exhibiting low levels of plastic deformation and viscoelastic behavior.2 This process is displayed in Figure 1 in the

form of a typical stress-strain curve for a material under uniaxial tension. The area under the linear unloading curve

is the recoverable, elastic strain energy density W and the remaining, shaded area under the curve is the dissipated

energy per unit initial volume WS.

Figure 1. Typical stress-strain curve displaying the elastic (W) and irrecoverable (WS ) potentials.

Both W and WS are functions of a set of internal state variables (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 yields the

thermodynamic driving force, fm, available for producing structural changes associated with the mth ISV.

fm =

∂Ws

∂Sm

(2)

It is shown in Ref. 24 that the total work potential is stationary with respect to each ISV.

∂WT

∂Sm

= 0 (3)

Additionally, Rice57 has shown that according to the second law of thermodynamics:

fmS ̇m ≥ 0 (4)

Equations (2), (3), and (4) form the foundation of a thermodynamically based work potential model for nonlinear

structural changes in a material. If one of the ISVs is chosen to be crack length, Equation 3 is analgous to Griffith’s

criterion commonly used in fracture mechanics.58

II.A. Dual-ISV formulation of ST to model progressive microdamage and transverse cracking

The generality of the evolution equations, presented in the previous section, allows for the response of a material

undergoing any number and type of structural changes to be captured. This is especially useful for modeling progressive damage in composites because multiple damage mechanisms arise during a typical loading history. For instance

3 of 22

NASA/TM—2012-217243 3
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.2, 3, 17, 54 The size of these flaws is typically on the order of that of the

fiber or smaller. Transverse cracks nucleate from preexisting flaws within the matrix but grow parallel to the fibers

and span the thickness of the lamina.5, 6, 8, 10, 18, 59, 60 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. Eventually, transverse cracking is succeeded by more catastrophic

damage mechanisms like interlamina delamination, fiber breakage, pull-out and bridging associated with macroscopic

laminate fracture.11, 47, 61

The average behavior of numerous coupon level laminates exhibiting progressive microdamage and transverse

cracking was successfully modeled in Refs. 2 and 54 using ST. A similar formulation will be used in this work, in

conjunction with the FEM, to model more complicated laminates. To incorporate microdamage and transverse cracking into ST, it is assumed that the inelastic work potential WS is a function of two ISVs: S which accounts for energy

dissipation via microdamage, and SC which captures available energy to form and advance transverse cracks. It can

be assumed that these ISVs are equal to the energy density required to advance the associated mechanism. Although

interaction between the mechanisms occurs through coupling in the evolution laws and adjustment of the stress and

strain fields with increasing damage, it is assumed that these ISVs are independent; in that, microdamage evolution

does not directly influence the growth of transverse cracks and vice versa. Thus, the total amount of dissipated energy

is a sum of the available energy used to progress each individual mechanism: WS = S + Sc. The total energy is then

given by

WT = W + S + SC (5)

Applying stationarity of potential energy (Equation (3)) to the definition of the total work potential yields two

evolution equations.

∂W

∂S = −1 (6)

∂W

∂SC

= −1 (7)

Additionally, invoking Equation (4) results in statements of the inadmissability of healing for both microdamage and

transverse cracking.

S ̇ ≥ 0 (8)

S ̇C ≥ 0 (9)

Equations (8) and (9) dictate that the amount of energy used to progress microdamage and transverse cracking can

never decrease; therefore, that energy has been dissipated into creating structural changes and cannot be recovered.

The combination of Equations (6)-(9) represent the evolution equations for microdamage and transverse cracking in

the matrix of the composite.

II.B. Application of dual-ISV ST to fiber reinforced laminae

II.B.1. 2-D plane stress constitutive law

To tailor the evolution equations presented in Section II.A to a unidirectional, fiber reinforced lamina, the stress-strain

relationships can be defined, under plane stress assumptions, in the principal material coordinates of the lamina as

σ11 = Q1111 + Q1222

σ22 = Q1211 + Q2222

τ12 = Q66γ12

(10)

where σij are the stress components, ij are the strain components, the 1-direction is aligned with the fibers in the

lamina, the ’2-direction is perpendicular to the fibers, γ12 is the engineering shear strain and Qij are the components

4 of 22

NASA/TM—2012-217243 4
of the stiffness matrix given by

Q11 =

E11

1 − ν12ν21

Q22 =

E22

1 − ν12ν21

Q12 = ν12Q22

Q66 = G12

ν21 =

ν12E22

E11

(11)

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.62 After assuming that the quantity ν12ν21 << 1,

Equations (11) become,

Q11 = E11

Q22 = E22

Q12 = ν12Q22

Q66 = G12

(12)

II.B.2. Determining the damage state

It is necessary to define the manner in which the moduli degrade as functions of the ISVs. Since the damage mechanisms associated with S and SC accrue only in the matrix of the composite, it is assumed that the moduli affected by

this damage are limited to E22 and G12. These moduli can be written as functions of S and SC .

E22 = E220es(S)ec(SC ) (13)

G12 = G120gs(S)gc(SC ) (14)

where E220 and G120 are the undamaged transverse and shear elastic moduli, es(S) and gs(S) are factors relating

the transverse and shear moduli to the microdamage S; ec(SC ) and gc(SC ) are functions relating the moduli to the

transverse cracking ISV SC . Sicking provided a procedure for determining the damage functions experimentally.2 The

experimental damage curves can then be fit with polynomials (such that values at S = 0 and SC = 0, are E22 = E220

and G12 = G120) and used in Equations (13) and (14).

The elastic strain energy density W

W =

Z

σijdij (15)

can be written using the plane stress constitutive relationships.

W =

1

2

(E11

2

11 + E22

2

22 + G12γ

2

12) + Q121122 (16)

Employing Equations (6) and (7) with (13), (14), and (16), while assuming the quantity Q12 = ν12Q22 (the effects

of evolving microdamage on the Poisson’s ratio have not be studied in detail) is constant (independent of S and SC ),

yields two ordinary differential equations which are be solved simultaneously for S and SC .



2

22

2

∂E22

∂S +

γ

2

12

2

∂G12

∂S = −1 (17)



2

22

2

∂E22

∂SC

+

γ

2

12

2

∂G12

∂SC

= −1 (18)

The above equations indicate that the microdamage and transverse cracking states depend only on the strain state,

elastic properties, and damage functions. Furthermore, Equations (17) and (18) preclude the use of any failure criteria

to determine the onset of matrix damage or failure. The damage functions, and slopes of the damage functions, dictate

the manner in which the ISVs evolve. Thus, for low strain levels, the damage ISVs remain small and have a negligible

effect on the response of the material. At intermediate strain levels the influence of the microdamage ISV becomes

more prevalent. At larger strain levels, the transverse cracking ISV begins to increase more rapidly. At even higher

strain levels, the ISVs are large enough that the transverse and shear stiffness begin to approach zero.

5 of 22

NASA/TM—2012-217243 5
Experimentally it has been determined that S and SC behave as 

3

, thus it is convenient to introduce reduced

damage variables Sr and Scr.

2

Sr ≡ S

1/3

(19)

Scr ≡ S

1/3

C

(20)

Furthermore, the use of reduced ISVs yields polynomial damage functions. With the reduced ISVs, the evolution

equations, Equations (17) and (18), become



2

22

2

E220

des (Sr)

dSr

ec (Scr) + γ

2

12

2

G120

dgs (Sr)

dSr

gc(Scr) = −3S

2

r

(21)



2

22

2

E220es (Sr)

dec (Scr)

dScr

+

γ

2

12

2

G120gs(Sr)

dgc (Scr)

dScr

= −3S

2

cr (22)

Once Sr and Scr are determined from Equations (21) and (22), while ensuring Equations (8) and (9), the transverse

and shear moduli can be degraded accordingly using Equations (13) and (14). The damage functions used in Equations

(13) and (14), and the procedures for obtaining, are given in Section III.

It is important to note that ST can account for compressive failure of a laminate resulting from fiber microbuckling

and kink band formation within the laminae. As a material point is loaded, the fibers begin to rotate. All calculations

are performed in the instantaneous material frame by tracking the fiber angle and transforming the field variables

accordingly. These rotations induce large, localized shear strains which, in turn degrade the matrix shear modulus,

allowing the fibers to rotate more easily. If the loading is compressive, this feedback will lead to a runaway instability

and regions of microbuckled, or kinked, fibers will develop. In previous studies,17 where the loading was compression

dominated, the instantaneous fiber rotation is captured in conjunction with ST. This has the advantage of predicting

fiber kinking failure, from a knowledge of basic constituent level properties, avoiding the use of an explicit fiber

direction compressive strength criterion, and the need for an additional coupon level test to measure compression

strength, Xc.

II.B.3. Decoupling the evolution equations

To obtain a solution for Equations (21) and (22) the damage functions must be explicitly defined. Typically fourthorder, or higher, polynomials are used for the microdamage and transverse cracking damage functions. The use of

these polynomial damage functions results in two coupled, bivariate polynomial evolution equations. The solution of

these equations produces numerous values for Sr and Scr. Only one set of solutions represents the physical behavior

of the material, the remaining solutions are mathematical artifacts. However to determine the correct solutions, all

solutions to the evolution equations must be obtained. This solution procedure can quickly become cumbersome and

computationally time consuming.

The solution of the evolution equations can be greatly simplified if the effects of microdamage and transverse

cracking on the transverse and shear stiffness are decoupled. Microdamage primarily manifests as shear microcracks

between fibers.3 A photograph of a [+45/−45]S composite specimen displaying these shear microcracks, taken under

a scanning electron microscope (SEM), is presented in Figure 2. Additionally, the influence of transverse cracking on

the shear modulus of a composite laminate is minimal.2 Thus, the evolution equations can be decoupled by assuming

that the shear modulus is affected only by microdamage, and the transverse stiffness is affected only by transverse

cracking.

E22 = E220ec(Scr) (23)

G12 = G120gs(Sr) (24)

Using these relationships, the decoupled evolution equations are

γ

2

12

2

G120

dgs(Sr)

dSr

= −3S

2

r

(25)



2

22

2

E220

dec(Scr)

dScr

= −3S

2

cr (26)

Even though the interaction between microdamage and transverse cracking is lost by decoupling the evolution equations, the solution procedure is greatly simplified making it more amenable for implementation into a finite element

6 of 22

NASA/TM—2012-217243 6
model where the evolution equations will need to be solved for every applicable material point at every time increment. As stated earlier, there is evidence to suggest that the uncoupled equations for damage evolution in laminated

composites is valid.2

Figure 2. SEM photograph showing shear microcracks in damaged [+45/ − 45]S specimen.

III. Experimental Program at NASA Langley Research Center

The response of IM7/8552 center-notched panels reinforced with buffer strip stiffeners was investigated at the

NASA Langley Research Center (LaRC) under the supersonics program.56 Each panel (Figure 3) consisted of a flat

12” x 34” composite laminate. The laminate was reinforced with two 1” wide buffer strips located 1.5” from each

edge. These buffer strips contained eight additional 0◦ plies not present in the rest of the laminate. Figure 4 displays

the cross section of the stiffened regions. The stacking sequences and thicknesses for the base laminate and buffer

strip regions are given in Table 1. The flat plate and buffer strips were co-cured as a single part. After the pristine

panel was cured, a centered notch that was 3” long and had notch tip diameter of 0.09735” was machined in the

panel. 5” x 12” steel grips were adhered to the ends panel. The grips were manufactured so they would conform to

the faces of the panel. One component of each set of grips was machined to accommodate the buffer strips, and the

other component was left flat. The bottom grip was clamped, preventing any displacements or rotations. The top grip

was secured to prevent any rotations or displacements in the x- or z-directions, while a tensile displacement of 5E-5

in./sec. was applied in the y-direction until the specimen failed catastrophically. Several strain gages, shown in Figure

3, were affixed to the specimen and were monitored during loading of the specimen along with the applied load and

edge displacement.

Region Stacking Thickness (in.)

Flat Panel [45◦/90◦/-45◦/0◦]S 0.044

Buffer Strip [45◦/90◦/-45◦/0o

5

]S 0.088

Table 1. Panel stacking sequence and thickness in flat and buffer strip-reinforced regions.

The buffer strips were added to the base laminate to increase the ultimate strength of the center-notched panels

without a substantial weight increase. Moreover, the intention of the buffer strips was to alter the failure mode of the

panel. Additional buffer strip configurations, not presented here, were also tested. These panels exhibited self similar

crack growth propagating from the notch tip towards the free edge. Upon reaching the buffer strips, however, the

direction of damage propagation turns so that the damage runs parallel to the buffer strip towards the gripped edge.

This behavior was not observed with the buffer strip configuration presented in this paper. Instead, the crack like

damage was arrested briefly at the buffer strip boundary before proceeding through the buffer strip towards the free

edge. Three specimens with the lay-up presented in Table 1 were tested in all.

The elastic properties for IM7/8552 are taken from Ref. 32 and given in Table 2. The shear microdamage curve was

obtained from a [+45/-45]5 tensile dogbone specimens. The normal transverse cracking damage curve was obtained

7 of 22

NASA/TM—2012-217243 7
Figure 3. Schematic of center-notched, buffer strip reinforced, tensile specimen.

Figure 4. Schematic of cross-section of buffer strip reinforced regions.

8 of 22

NASA/TM—2012-217243 8
from multi-angle AS4/3502 specimens, presented in Ref. 2. The coefficients of this curve were scaled by the ratio of

the elastic E220 for IM7/8552 to the elastic E220 for AS4/8502.

e

IM7/8552

ci =

EIM7/8552

220

EAS4/3502

220

e

AS4/3502

ci , i 6= 0 (27)

This scaling procedure was used in the successful failure prediction of notched compression specimens as reported

in Ref. 63. Figures 5 and 6 show the shear microdamage curve and normal transverse cracking damage curve as

functions of the respective reduced ISVs. The coefficients of the polynomial curves are presented in Table 3.

Property Value

E11 (msi) 24.86

E22 (msi) 1.78

ν12 0.32

G12 (msi) 0.68

Table 2. Elastic properties and axial strain allowable for IM7/8552.

Figure 5. Shear microdamage function.

Figure 6. Normal transverse cracking damage function.

9 of 22

NASA/TM—2012-217243 9
Function c0 c1 c2 c3 c4 c5

gs(Sr) 1.0000 -6.3955E-2 -4.7716E-2 1.0566E-2 -8.3609E-4 2.2997E-5

ec(Scr) 1.0000 4.9809E-2 -6.7644E-2 -1.5854E-2 2.2461E-3 0.0000

Table 3. Polynomial coefficients of damage curves.

IV. Finite Element Model

IV.A. Finite element discretization and boundary conditions

One quarter of the center-notched, buffer-strip panel, presented in Section III, was discretized using ABAQUS, S4R,

4-noded, quadrilateral and S3R, 3-noded, triangular, two-dimensional (2D), plane stress, shell elements.64 This mesh

contained 9257 S4R and 42 S3R elements, with a high density of elements near the notch tips to the free edge of the

panel. The high density regions were chosen to accurately capture the fields near the notch tip and along the expected

axial failure path. This mesh is displayed in Figure 7. The portions highlighted in green indicate buffer strip-reinforced

regions.

Figure 7. Finite element mesh of quarter of buffer strip panel domain. Buffer strip reinforced elements highlighted in red.

Displacements in the x-, y-, and z-directions as well as all rotations along the bottom edge of all the panel are

constrained. Similarly for the top edges, the displacements in the x- and z-directions and all rotations are prohibited.

A uniform displacement in the y-direction is applied dynamically at a rate of 0.886 in./sec. up to 0.0886 in. over

a duration of 0.1 sec. Initially, mass scaling is introduced into the model on an element by element basis, using the

ABAQUS/Explicit command *MASS SCALING, to increase the minimum stable time increment of the problem to

1.0E-5 sec. The input deck card *VARIABLE MASS SCALING is also used to scale the mass of the elements every

1000 time increments if necessary.

IV.B. In-plane damage and failure

Dual ISV ST is implemented into the FEM model using the commercially available ABAQUS 6.8-1 finite element software.64 The ST constitutive laws are integrated into ABAQUS/Explicit through a user material subroutine (VUMAT).

10 of 22

NASA/TM—2012-217243 10
Axial failure is also included in the VUMAT to capture the effects of fiber breakage and other axial failure mechanisms.

Prior to calculating the damage state of the material point, an axial maximum strain failure criterion is evaluated,

ef =

11



U

11

(28)

where 11 is the strain in the fiber direction of the layer and 

U

11 is the ultimate strain. If ef ≥ 0 the element is deleted

and removed from subsequent calculations. If axial failure has not occurred; then, at each integration point in an

element obeying ST, the local strain fields are used in Equations (25) and (26) to calculate Sr and Scr. Equations (8)

and (9) are enforced manually; if an ISV decreases, the value of the ISV at the end of the previous increment is used.

Once the damage state is determined, the transverse and shear stiffness are updated using Equations (23) and (24).

The stress field is then calculated using the degraded moduli and Equation (10).

The IM7/8552 lamina elastic properties are given in Table 2, and the damage function coefficients are presented

in Table 3. It should be mentioned that the axial failure criterion, Equation (28) and the degradation scheme used are

very mesh dependent. The experimental, average, maximum value of strain near the notch (Sg-8) was used to estimate

the maximum strain allowable for the simulation. This strain allowable was 2992 microstrain (μ). Due to the very

steep strain gradients near the notch, the average strain in a region the size of the element used at the notch tip in the

simulation would be far less than the strain reported by the strain gage. Therefore, the experimentally obtained strain

allowable was reduced by 25% to arrive at an average allowable for the axial failure strain of the elements used in

the model: 

U

11 = 2250μ. Ongoing studies will include a sensitivity analysis of 

U

11 to examine how it influences the

overall results of the model. Additionally, future models will incorporate more sophisticated axial failure models that

can account for pullout and friction during axial failure of a fiber composite,65 as well as possible non-local effects,66

thus, alleviating the dependence of axial failure on mesh size and density.

V. Results

V.A. Quantitative results compared to experimental data

Three center-notched, buffer strip-reinforced panels were tested under uniaxial tension at NASA LaRC. Load versus

edge displacement as well as local strain gage data was recorded. The location of the strain gages is displayed in

Figure 3. Each strain gage shown in the diagram represents multiple strain gages placed in equivalent locations on the

panel. The data presented here is the average value of those analogous strain gage readings.

A comparison of the global load-deflection curve obtained from the model to those from experiment is plotted in

Figure 8. The simulation is able to capture both the linear and non-linear load versus displacement response of the

panel. There was significant variability between experimental results for the nominally same panel. The ultimate loads

(P

ult) for the three panels were 37,407 lb., 30,451 lb., and 23,471 lb. Thus, the average ultimate load for the three

tests was 30,443 lb. with a standard deviation of 6,968 lb. The finite element model yielded an ultimate load of 32,844

which is approximately 8% greater than the average. This indicates that the 25% reduction imposed on the average

maximum strain near the notch tip to obtain the axial failure strain allowable was slightly non-conservative, or that

energy was dissipated via other damage mechanisms not included in the finite element model, such as delamination or

fiber bridging. Discrepancies between the predicted ultimate load and the experimental values could also be a result

of variability in the experiment, or a artifact of the degradation scheme used subsequent to satisfaction of the failure

criterion. The ultimate loads are summarized in Table 4. In addition to an adequate ultimate load prediction, the

simulation also exhibits several discontinuous load-drops, that are seen in the experimental results, prior to reaching

the ultimate load. This may indicate that the failure and damage mechanism evolution predicted with the finite element

model are accurate. An overview of these damage and failure mechanisms are presented in the next section.

Figure 9 shows the ratio of the kinetic energy (KE) to the total energy (ET otal) of the system as a function of time.

The kinetic energy remains less than 5% of the total energy throughout the simulation except at the very beginning of

the model (well within the linear regime), and just prior to ultimate failure. This indicates that the mass scaling and

rapid loading rate used to speed up the total simulation time are not influencing the results in a non-physical manner.

Local strain gage data for the nine gages (see Figure 3) are plotted in Figure 10. The linear stiffness of the numerical

model is not in correspondence with any of the strain gage stiffness’. It should be noted that there is a discrepancy

between the stiffness of the experimental load-displacement curve and the load-axial strain curve of Sg-1 far away

from the notch. The reason for this discrepancy has yet to be uncovered; however, the trends predicted by the finite

element model can still be compared to those observed in the experiment. Since the ultimate load obtained from the

model was closest to that in Test2, local model data will be compared most closely to that test.

11 of 22

NASA/TM—2012-217243 11
Figure 8. Experimental and numerical load vs. displacement data for buffer strip panel.

Figure 9. Kinetic energy history.

12 of 22

NASA/TM—2012-217243 12
(a) Sg-1 (b) Sg-2 (c) Sg-3

(d) Sg-4 (e) Sg-5 (f) Sg-6

(g) Sg-7 (h) Sg-8 (i) Sg-9

Figure 10. Experimental and numerical load vs. axial strain data for buffer strip panel.

13 of 22

NASA/TM—2012-217243 13
At Sg-1, the initial load at which the strain begins to relax in the numerical simulation is comparable to the load in

Test2. This strain relaxation far from the notch is a result of damage occuring near the notch tip. The strain response

at Sg-2, which is located on the buffer strip far from the notch, is linear in both the simulation and the experiment.

Sg-3, which is located at the boundary of the buffer strip but 6 in. above the notch, exhibits an initial strain relaxation

in both Test2 and the finite element model at the same load level as observed in Sg-1. Jumps in the strain levels at Sg-4

and Sg-6 occur in the model at global loads significantly less than in Test2. This indicates that the damage and failure

evolution progresses through the buffer strip more rapidly in the numerical model than in the experiment. However,

data for Sg-6, which lies on the buffer strip directly in front of the notch, matches well with the experiment. The

same can be seen at Sg-7. Experimental data at the notch tip (Sg-8) displays a high level of non-linearity. This is

not observed in the FE model and may be a facet of the extremely steep strain gradients at the notch tip which would

require an extremely fine mesh at that region. Similarly, the numerical model does not capture the non-linearity shown

in the experiment at Sg-9 in the center of the panel directly above the notch. A relatively large element is present in

the model at that location and is used to represent that strain gage; a more refined mesh in that region may capture the

appropriate response, and this is currently being studied. Additionally, poor performance of the FEM model at Sg-8

and -9 may be due to the lack of transverse-shear damage coupling.

V.B. Analysis damage and failure mechanisms observed in numerical simulation

The results of the finite element model indicate that large degrees of microdamage, transverse cracking, and axial

failure are confined to regions in front of the notch tip. Thus, Sg-1, Sg-2, and Sg-3 (see Figure 3), which lie in regions

far from the notch, should not be directly influenced by damage or failure, but rather are affected through redistribution

of the strain fields. Figure 11 shows the evolution of the microdamage and transverse crack ISVs with respect to the

applied displacement, in each of the top four plies, at the strain gage locations near and in front of the notch tip (Sg4-Sg-9). The ISVs are normalized by their maximum values (S

max

r = 12.08 psi1/3

and S

max

cr = 3.52 psi1/3

) so that

a value of one in Figure 11 corresponds to complete degradation of the appropriate modulus (see Equations (23) and

(24)).

Similar damage evolution can be seen at Sg-4, Sg-5, and Sg-6. Microdamage evolves progressively, and equally,

in the +45◦and -45o plies. Very low levels of Sr arise at these gages in the 0◦ plies, and no microdamage is observed

in the 90◦ ply until catastrophic failure. Transverse crack evolution remains dormant in all plies, until after a 0.13 in.

edge displacement is achieved; after which the degree of damage due to transverse cracking rises rapidly in the 90◦

ply until catastrophic failure. This increase in Scr is associated with propagation of axial failure from the notch tip,

as is evidenced in the proceeding damage contour analysis. In the +45◦, -45◦, and 0◦ plies, Scr remains zero until

catastrophic failure.

At the strain gage on the boundary of the buffer strip closest to the notch, only very low levels of microdamage

are observed in the +45◦and -45◦ plies. However, both Sr and Scr jump suddenly at an edge displacement of 0.113”.

This jump in the ISVs is associated with axial failure evolving from the notch meeting the buffer strip boundary.

At the notch tip (Sg-8), transverse cracking in the +45◦and 90◦layers develop rapidly at low displacement levels.

Expeditious advancement in Scr is accompanied by a steady increase of Sr in the +45◦and -45◦layers near the notch

tip. Less severe microdamage evolution also occurs in the 90◦and 0◦laminae. At a global displacement of 0.1”. there

is a large jump in Scr in the -45◦layer. Scr remains zero in the 0◦ plies; this is most likely a result of axial failure of

these layers at the notch tip before transverse cracks can develop. Above the notch at Sg-9, low levels of microdamage

in the +45◦and -45◦are observed. Overall, the evolution of microdamage in the plies is progressive; whereas, the

advancement of transverse cracks is sudden and discontinuous.

Figures 12-15 show matrix microdamage and transverse cracking contours, as well as axial failure paths, at several

Load (lb.) Deviation from

Experimental Average

Test1 37,407 23%

Test2 30,451 0.0003%

Test3 23,471 23%

Average 30,443 N/A

FEM 32,884 8%

Table 4. Ultimate load information for experiments and finite element model.

14 of 22

NASA/TM—2012-217243 14
(a) Sg-4 (b) Sg-5 (c) Sg-6

(d) Sg-7 (e) Sg-8 (f) Sg-9

Figure 11. Microdamage and transverse crack evolution at strain gage locations in front notch tip.

15 of 22

NASA/TM—2012-217243 15
(a) Sr, t = 0.0046s (b) Sr, t = 0.0049334s (c) Sr, t = 0.005334s

(d) Scr, t = 0.0046s (e) Scr, t = 0.0049334s (f) Scr, t = 0.005334s

(g) Sr, t = 0.0054667s (h) Sr, t = 0.0061334s (i) Sr, t = 0.0064001s

(j) Scr, t = 0.0054667s (k) Scr, t = 0.0061334s (l) Scr, t = 0.0064001s

Figure 12. Microdamage, transverse crack evolution and axial failure in top 45◦ply.

16 of 22

NASA/TM—2012-217243 16
(a) Sr, t = 0.0046s (b) Sr, t = 0.0049334s (c) Sr, t = 0.005334s

(d) Scr, t = 0.0046s (e) Scr, t = 0.0049334s (f) Scr, t = 0.005334s

(g) Sr, t = 0.0054667s (h) Sr, t = 0.0061334s (i) Sr, t = 0.0064001s

(j) Scr, t = 0.0054667s (k) Scr, t = 0.0061334s (l) Scr, t = 0.0064001s

Figure 13. Microdamage, transverse crack evolution and axial failure in top 90◦ply.

17 of 22

NASA/TM—2012-217243 17
(a) Sr, t = 0.0046s (b) Sr, t = 0.0049334s (c) Sr, t = 0.005334s

(d) Scr, t = 0.0046s (e) Scr, t = 0.0049334s (f) Scr, t = 0.005334s

(g) Sr, t = 0.0054667s (h) Sr, t = 0.0061334s (i) Sr, t = 0.0064001s

(j) Scr, t = 0.0054667s (k) Scr, t = 0.0061334s (l) Scr, t = 0.0064001s

Figure 14. Microdamage, transverse crack evolution and axial failure in top -45◦ply.

18 of 22

NASA/TM—2012-217243 18
(a) Sr, t = 0.0046s (b) Sr, t = 0.0049334s (c) Sr, t = 0.005334s

(d) Scr, t = 0.0046s (e) Scr, t = 0.0049334s (f) Scr, t = 0.005334s

(g) Sr, t = 0.0054667s (h) Sr, t = 0.0061334s (i) Sr, t = 0.0064001s

(j) Scr, t = 0.0054667s (k) Scr, t = 0.0061334s (l) Scr, t = 0.0064001s

Figure 15. Microdamage, transverse crack evolution and axial failure in 0◦plies.

19 of 22

NASA/TM—2012-217243 19
key time points in the top four layers (45◦

, -45◦, 90◦, and 0◦) of the panel acreage and buffer strip. These snapshots

depict the evolution of, as well as the interactions between, the three damage mechanisms modeled in this work. The

damage variables, Sr and Scr are represented continuously in the colored contours. Blue indicates zero, or very little,

damage, while red indicates the maximum damage (i.e. when the stiffness has degraded to zero). Turquoise, green,

yellow, and orange bridge the scale between minimum damage (blue) and maximum damage (red). Elements that have

satisfied the axial failure criterion, or exceed some excessively large strain value, are considered to not contribute any

stiffness to the panel and are removed. Elements in the buffer strip region are outlined in green to display the effects

of the pad-up in this region on the crack path and damage evolution.

In the outer 45◦layer, prior to initiation of axial failure, a very small amount of microdamage is present near the

notch tip and is accompanied by a localized region of maximum transverse cracking damage. Shortly after this, axial

failure occurs in these elements. When the maximum strain criterion is satisfied in an element, that element is deleted,

and the corresponding void that develops essentially represents a macroscopic crack in the material. Additionally, any

elements that have reached S

max

cr , and have a transverse strain that has surpassed the axial strain allowable, are also

removed to prevent excessive element distortion in the simulation. This is a reasonable assumption since transverse

failure strains are typically lower than axial failure strains, and an element that has reached S

max

cr has most likely been

loaded predominantly in the transverse direction and is not contributing axially. When the crack reaches the buffer

strips it changes direction and runs parallel to the buffer strip. The crack is then arrested for a brief moment before

it advances laterally into the buffer strip (Figure 12(c)). Once inside the buffer strip, the crack turns again and runs

along the buffer strip. The subsequent point of arrest for the crack is in a transition region between a very dense mesh

and a coarser mesh. It may be that this transition region is influencing the behavior of the crack in a non-physical

way; however, similar behavior was observed in regions with a uniformly dense mesh (see Figures 12(c) and 12(f)).

A crack branches from the original crack and propagates perpendicular to the length of the buffer strip. Shortly after

this branched crack forms, numerous other branched cracks advance through the buffer strip leading to catastrophic

failure of the panel. The evolution of the crack in the 45◦layer is not accompanied by much microdamage. However,

large regions exhibiting high levels of transverse cracking surround the macroscopic crack is evident.

The same macroscopic crack progression is observed in the top 90◦, -45◦, and 0◦layers. Additionally, the accumulation in all four plies is very similar. The 0◦

and 90◦ plies exhibit more microdamage encompassing a larger area

surrounding the crack as it advances. The 90◦layer exhibits widespread damage due to transverse cracking (Figure

13. In the 0◦layer, high values of Scr are limited to the vicinity surrounding the axial failure path.

The individual damage mechanisms (i.e. microdamage, transverse cracking, and axial failure) are seen to be

interactive. Accumulation of Sr and Scr near the notch tip leads to a local transfer of load into the 0◦ ply from

adjacent layers, which subsequently leads to local axial failure. Similarly, the fields surrounding the advancing crack

lead to rapid matrix degradation and further crack advancement. Upon, reaching the buffer strip, the path of least

resistance to axial failure lies along the longitudinal direction of the buffer strip. However, substantial regions of

severe matrix degradation in the non-0◦ plies in the buffer strip provide an avenue for axial failure to occur within

the 0◦layers of the buffer strip through load transfer. Once the reinforcement provided by the buffer strip begins to

diminish, catastrophic failure is imminent and this is captured in the present work.

VI. Conclusions

A finite element based implementation of progressive damage accumulation and failure, or in short, progressive

failure analysis (PFA), that is capable of capturing three predominant damage mechanisms in an FRL (namely: matrix

shear microdamage, matrix transverse cracking, and axial failure) has been presented. The matrix damage mechanisms

are modeled using a dual ISV formulation of ST, which accounts for the damage development mechanisms via the

energy required to create structural changes associated with the damage modes in the material. A maximum strain

cut-off is chosen to capture axial failure within layers of a laminate. The capabilities of the model are demonstrated

by analyzing a center-notched, buffer-strip reinforced composite panel, for which there is experimental data. Good

quantitative agreement with global load-displacement, and local strain gage data is exhibited. Moreover, the resulting damage evolution in the model provides insight into the interactions between the individual damage mechanisms

within the composite that lead to the global observed failure modes. To obtain subtle details, transparency, and accessibility into the evolution of damage modes in a complex composite structure with experimentation alone would be

difficult.

While a significant number of important damage accumulation and failure mechanisms have been captured, delamination failure has not been included in the model, but it can be incorporated by including DCZM elements at the

interfaces between adjacent shell or three-dimensional sublaminate layers. Further investigations will focus on the

20 of 22

NASA/TM—2012-217243 20
influence of interlamina damage on the response of the buffer strip panel. Additionally, the axial failure theory used

does not properly capture the physics involved in axial damage of FRLs. A sensitivity analysis will be performed to

determine the influence of the axial failure strain on the ultimate load of, and failure path within, the panel. Future

efforts will also attempt to forgo the maximum strain criterion for axial failure in favor of a physically-based, non-local

axial damage theory,66 or a complete micromechanics simulation of axial failure that includes pull-out and frictional

slip.65

References

1Beaumont, P. W. R., Dimant, R. A., and Shercliff, H. R., “Failure processes in composite materials: getting physical,” J. Mater. Sci., Vol. 41,

2006, pp. 6526–6546.

2Sicking, D. L., Mechanical Characterization of Nonlinear Laminated Composites with Transverse Crack Growth, Ph.D. thesis, Texas A&M

University, College Station, TX, 1992.

3Ng, W. H., Salvi, A. G., and Waas, A. M., “Characterization of the in-situ non-linear shear response of laminated fiber-reinforced composites,” Compos. Sci. Technol., Accepted February 24, 2010.

4Dvorak, G. J. and Laws, N., “Analysis of progressive matrix cracking in composite laminates I. Thermoelastic properties of a ply with

cracks,” J. Compos. Mater., Vol. 19, 1985, pp. 216–234.

5Talreja, R., “Transverse cracking and stiffness reduction in composite laminates,” J. Composite Materials, Vol. 19, 1985, pp. 355–275.

6Allen, D. H., Harris, C. E., and Groves, S. E., “A thermomechanical constitutive theory for elastic composites with distributed damage - II.

Aplication to matrix cracking in laminated composites,” Int. J. Solids Struct., Vol. 23, No. 9, 1987, pp. 1319–1338.

7Tan, S. C. and Nuismer, R. J., “A theory for progressive matrix cracking in composite laminates,” J. Compos. Mater., Vol. 23, 1989,

pp. 1029–1047.

8Noda, J., Okabe, T., Takeda, N., and Shimizu, M., “Tensile strength of CFRP cross-ply laminates containing transverse cracks,” Adv. Comps.

Mater., Vol. 15, No. 1, 2006, pp. 81–93.

9Varna, J., Joffe, R., and Talreja, R., “A synergistic damage-mechanics analysis of transverse cracking in [±θ/904]S

,” Compos. Sci. Technol.,

Vol. 61, 2001, pp. 657–665.

10Yang, Q. and Cox, B., “Cohesive models for damage evolution in laminated composites,” Int. J. Frac., Vol. 133, 2005, pp. 107–137.

11McCartney, L. N., “Mechanics for the growth of bridged cracks in composite materials: Part I. Basic Principles,” J. Compos. Technol. &

Res., Vol. 14, No. 3, 1992, pp. 133–154.

12Rudraraju, S. S., Salvi, A., Graikipati, K., and Waas, A. M., “In-plnae fracture of laminated fiber reinforced composites with varying fracture

resistance: experimental obesrvations and numerical crack propagation simulations,” 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural

Dynamics, and Materials Conference, 4-7 May 2009.

13Budiansky, B. and Fleck, N. A., “Compressive failure of fiber composites,” J. Mech. Phys. Solids, Vol. 41, No. 1, 1993, pp. 183–211.

14Fleck, N. A., Deng, L., and Budiansky, B., “Prediction of kink width in compressed fiber composites,” J. App. Mech., Vol. 62, 1995,

pp. 329–337.

15Schapery, R. A., “Prediction of compressive strength and kink bands in composites using a work potential,” Int. J. Solids Structures, Vol. 32,

No. 6, 1995, pp. 739–765.

16Yerramalli, C. S. and Waas, A. M., “A failure criterion for fiber reinforced polymer composites under combined compression-torsion loading,” Int. J. Solids Structures, Vol. 40, No. 5, 2003, pp. 1139–1164.

17Basu, S., Waas, A. M., and Ambur, D. R., “Compressive failure of fiber composites under multiaxial loading,” Int. J. Solids Structures,

Vol. 44, No. 9, 2006, pp. 2648–2676.

18Green, B. G., Wisnom, M. R., and Hallett, S. R., “An experimental investigation into the tensile strength scaling of notched composites,”

Compos.: Part A, Vol. 38, 2007, pp. 867–878.

19Pierron, F., Green, B., and Winsom, M. R., “Full-field assessment of the damage process of laminated composite open-hole tensile specimens.

Part I: Methodology,” Composites: Part A, Vol. 38, 2007, pp. 2307–2320.

20Pierron, F., Green, B., and Winsom, M. R., “Full-field assessment of the damage process of laminated composite open-hole tensile specimens.

Part II: Experimental results,” Composites: Part A, Vol. 38, 2007, pp. 2321–2332.

21Waas, A. M., C. D. Babcock, J., and Knauss, W. G., “An Experimental Study of Compression Failure of Fibrous Laminates in the Presence

of Stress Gradients,” Int. J. Solids Struct., Vol. 26, No. 9, 1990.

22Khamesh, A. and Waas, A. M., “Failure of uniply model composites under compression,” ASME Trans. J. of Eng. Mater. Technol., Vol. 114,

1992, pp. 304–310.

23Waas, A. M., C. D. Babcock, J., and Knauss, W. G., “An Experimental Study of Compression Failure of Fibrous Laminates in the Presence

of Stress Gradients,” Int. J. Solids Struct., Vol. 26, No. 9, 1990.

24Schapery, R. A., “A theory of mechanical behaviour of elastic media with growing damage and other changes in structure,” J. Mech. Phys.

Solids, Vol. 38, No. 2, 1990, pp. 1725–1797.

25Lemaitre, J. and Chaboche, J.-L., Mechanics of Solid Materials, Cambridge University Press, 1994.

26Talreja, R., “Matrix microcracking,” Damage characterization by internal variables, edited by R. Talreja, Elsevier Science B.V., 1994.

27Ladeveze, P., “Inelastic strains and damage,” Damage Mechanics of Composite Materials, edited by R. Talreja, Elsevier Science B.V., 1994.

28Lemaitre, J., A Course on Damage Mechanics, 2nd Ed., Springer-Verlag, Berlin, Heidelberg, New York, 1996.

29Sukumar, N., Moes, N., Moran, B., and Belytschko, T., “An extended finite element method (X-FEM) for two- and three-dimensional crack  ̈

modeling,” Int. J. Numer. Methods Eng., Vol. 48, No. 11, 2000, pp. 1741–1760.

30McCartney, L. N., “Physically based damage models for laminated composites,” Proc. Instn Mech. Engrs, Vol. 217, IMechE, 2003, pp.

163–199.

21 of 22

NASA/TM—2012-217243 21
31Li, S., Thouless, M. D., Waas, A. M., Schroeder, J. A., and Zavattierid, P. D., “Use of a cohesive zone model to analyze the frcture of a

fber-reinfoced polymermatrix composite,” Composites Science and Technology, Vol. 65, 2005, pp. 537–549.

32Camanho, P. P., Maimi, P., and Davila, C. G., “Prediction of size effects in notched laminates using continuum damage mechanics,” Compos.

Sci. and Technol., Vol. 67, 2007, pp. 2715–2727.

33Schuecker, C. and Pettermann, H. E., “Fiber reinforced laminates: progressive damage modeling based on failure mechanisms,” Arch.

Comput. Methods Eng., Vol. 15, 2008, pp. 163–184.

34Aboudi, J., Mechanics of Composite Materials: A Unified Micromechanical Approach, Elsevier Amsterdam, 1991.

35Paley, M. and Aboudi, J., “Micromechanical analysis of composites by the generalized cells model,” Mechanics of Materials, Vol. 14, 1992,

pp. 127–139.

36Bazant, Z. P., “Nonlocal damage theory based on micromechanics of crack interactions,” J. Eng. Mech. - ASCE.

37Prat, P. C. and Bazant, Z. P., “Tangential stiffness of elastic materials with systems of growing or closing cracks,” J. Mech. Phys. Solids,

Vol. 45, No. 4, 1997, pp. 611–636.

38Totry, E., Gonzalez, C., and Llorca, J., “Influence of the loading path on the strength of fiber-reinforced composites subjected to transverse

compression and shear,” Int. J. Solids Structures, Vol. 45, 2008, pp. 1663–1675.

39Bednarcyk, B. A., Aboudi, J., and Arnold, S. M., “Micromechanics modeling of composite materials subjected to multiaxial progressive

continuum damage,” 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 4-7 May 2009.

40Kouznersova, V., Brekelmans, W. A. M., and Baaijens, F. P. T., “An approach to micro-maro modeling of heterogenerous materials,” Comput.

Mech., Vol. 27, 2001, pp. 37–48.

41Bednarcyk, B. A. and Arnold, S. M., “A framework for performing multiscale stochastic progressive failure analysis of composite structures,”

Proceedings of the 2006 ABAQUS User’s Conference, 23-25, May 2006.

42Talreja, R., “Multi-scale modeling in damage mechanics of composite materials,” J. Mater. Sci., Vol. 41, 2006, pp. 6800–6812.

43Kwon, Y. W., Allen, D. H., and Talreja, R., editors, Multiscale Modeling and Simulation of Composite Materials and Structures, Springer,

2008.

44Talreja, R., “Multi-scale modeling of composite solids with damage,” Proceedings of the American Society for Composites Twenty-third

Technical Conference, 9-11, Sept. 2008.

45Pineda, E. J., Waas, A. M., Bednarcyk, B. A., and Collier, C. S., “A novel, multiscale high fidelity progressive damage and failure modeling

approach for laminated fiber reinforced composites,” Proceedings of the American Society for Composites Twenty-third Technical Conference, 7-10,

Apr. 2008.

46Turon, A., Camanho, P. P., Costa, J., and Davila, C. G., “A damage model for the simulation of delamination in advanced composites under

variable-mode loading,” Mech. Mater., Vol. 38, No. 11, 2006.

47Hallett, S. R., Jiang, W.-G., Khan, B., and Winsom, M. R., “Modelling the interaction between matrix cracks and delamination damage in

scaled quasi-isotropic specimens,” Compos. Sci. and Technol., Vol. 68, 2008, pp. 80–89.

48Guiamatsia, I., Anderson, J. K., Davies, G. A. O., and Iannucci, L., “Decohesion finite element with enriched basis functions for delamination,” Compos. Sci. Technol., Vol. 69, 2009, pp. 2616–2624.

49Song, S. J. and Waas, A. M., “An energy based model for mixed mode failure of laminated composites,” AIAA Journal, Vol. 33, No. 4, 1995.

50Xie, D. E., Salvi, A., Sun, C. E., Waas, A. M., and Caliskan, A., “Discrete cohesive zone model to simulate static fracture in 2-D triaxially

braided carbon fiber composites,” J. Compos. Mater., Vol. 40, No. 22, 2006.

51Xie, D. and Waas, A. M., “Discrete cohesive zone model for mixed-mode fracture using finite element analysis,” Engineering Fracture

Mechanics, Vol. 73, No. 13, 2006, pp. 1783–1796.

52Gustafson, P. A. and Waas, A. M., “The influence of adhesive constitutive parameters in cohesive zone finite element models of adhesively

bonded joints,” Int. J. Solids Struct., 209, pp. 2201–2215.

53Ranatunga, V., Bednarcyk, B. A., and Arnold, S. M., “Modeling delamination in composites via continuum interfacial displacement discontinuities,” 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 4-7 May 2009.

54Schapery, R. A. and Sicking, D. L., “A theory of mechanical behaviour of elastic media with growing damage and other changes in structure,”

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

55Pineda, E. J., Waas, A. M., Bednarcyk, B. A., Collier, C. S., and Yarrington, P. W., “Progressive damage and faiure modeling in notched

laminated fiber reinforced composites,” Int. J. Fract., Vol. 158, 2009, pp. 125–143.

56Bogert, P. B. and Satyanaranya, A., private communication, 2010.

57Rice, J. R., “Inelastic constitutive relations for solids: an internal-variable theory and its application to metal plasticity,” J. Mech. Phys.

Solids, Vol. 19, 1971, pp. 433–455.

58Sanford, R. J., Principles of Fracture Mechanics, Pearson Education, Inc., 2002.

59Laws, N. and Dvorak, G. J., “Progressive transverse cracking in composite laminates,” J. Compos. Mater., Vol. 22, 1988, pp. 900–916.

60Gudmundson, P. and Ostlund, S., “Numerical verification of a procedure for calculation of elastic constants in microcracking composite

laminates,” J. Compos. Mater., Vol. 26, 1992, pp. 2480–2492.

61McCartney, L. N., “Theory of stress transfer in a 0o

-90o-0

o

cross-ply laminate containing a parallel array of transverse cracks,” J. Mech.

Phys. Solids, Vol. 40, No. 1, 1992, pp. 27–68.

62Jones, R. M., Mechanics of Composite Materials, 2nd Ed., Taylor and Francis, Inc., Philadelphia, 1999.

63Basu, S., Computational Modeling of Progrssive Failure and Damage in Composite Laminates, Ph.D. thesis, University of Michigan, Ann

Arbor, MI, 2005.

64SIMULIA, Providence, RI, ABAQUS User’s Manual, Vol. 1-3, Version 6.8-1, 2008.

65Rudraraju, S., Salvi, A., Garikipati, K., and Waas, A. M., “Micromechanics based full numerical simulation of axial failure in fiber composites,” in preparation, 2010.

66Pineda, E. J., Waas, A. M., Bednarcyk, B. A. ., and Arnold, S. M., “Non-local effects in axial failure of fiber composites,” in preparation,

2010.

22 of 22

NASA/TM—2012-217243 22

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-05-2012 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Computational Implementation of a Thermodynamically Based Work Potential Model For 

Progressive Microdamage and Transverse Cracking in Fiber-Reinforced Laminates 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Pineda, Evan, J.; Waas, Anthony, M.; Bednarcyk, Brett, A.; Collier, Craig, S. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 727950.04.04.22 

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-17999 

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-2012-217243 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Category: 39 

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 continuum-level, dual internal state variable, thermodynamically based, work potential model, Schapery Theory, is used capture the 

effects of two matrix damage mechanisms in a fiber-reinforced laminated composite: microdamage and transverse cracking. Matrix 

microdamage accrues primarily in the form of shear microcracks between the fibers of the composite. Whereas, larger transverse matrix 

cracks typically span the thickness of a lamina and run parallel to the fibers. Schapery Theory uses the energy potential required to advance 

structural changes, associated with the damage mechanisms, to govern damage growth through a set of internal state variables. These state 

variables are used to quantify the stiffness degradation resulting from damage growth. The transverse and shear stiffness’ of the lamina are 

related to the internal state variables through a set of measurable damage functions. Additionally, the damage variables for a given strain 

state can be calculated from a set of evolution equations. These evolution equations and damage functions are implemented into the finite 

element method and used to govern the constitutive response of the material points in the model. Additionally, an axial failure criterion is 

included in the model. The response of a center-notched, buffer strip-stiffened panel subjected to uniaxial tension is investigated and results 

are compared to experiment. 

15. SUBJECT TERMS

Composites; Progressive failure analysis; Continuum damage; Laminates 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

30 

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