Evan J. Pineda and Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

Steven M. Arnold

Glenn Research Center, Cleveland, Ohio

Progressive Failure of a Unidirectional

Fiber-Reinforced Composite Using the

Method of Cells: Discretization Objective

Computational Results

NASA/TM—2012-217649

August 2012

AIAA–2012–1539
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c 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 Offi cer. 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 signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c 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. Scientifi c 

and technical fi ndings 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. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c 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 to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443–757–5803

• Phone the NASA STI Information Desk at

 443–757–5802

• Write to:

 STI Information Desk

 NASA Center for AeroSpace Information

 7115 Standard Drive

 Hanover, MD 21076–1320
Evan J. Pineda and Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

Steven M. Arnold

Glenn Research Center, Cleveland, Ohio

Progressive Failure of a Unidirectional

Fiber-Reinforced Composite Using the

Method of Cells: Discretization Objective

Computational Results

NASA/TM—2012-217649

August 2012

AIAA–2012–1539

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

53rd Structures, Structural Dynamics, and Materials (SDM) Conference

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

Honolulu, Hawaii, April 23–26, 2012
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 identifi cation 

only. Their usage does not constitute an offi cial endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

Level of Review: This material has been technically reviewed by technical management. 
The smeared crack band theory is implemented within the generalized method of cells

and high-fidelity generalized method of cells micromechanics models to capture progressive

failure within the constituents of a composite material while retaining objectivity with

respect to the size of the discretization elements used in the model. An repeating unit

cell containing 13 randomly arranged fibers is modeled and subjected to a combination

of transverse tension/compression and transverse shear loading. The implementation is

verified against experimental data (where available), and an equivalent finite element model

utilizing the same implementation of the crack band theory. To evaluate the performance

of the crack band theory within a repeating unit cell that is more amenable to a multiscale

implementation, a single fiber is modeled with generalized method of cells and high-fidelity

generalized method of cells using a relatively coarse subcell mesh which is subjected to the

same loading scenarios as the multiple fiber repeating unit cell. The generalized method

of cells and high-fidelity generalized method of cells models are validated against a very

refined finite element model.

I. Introduction

Micromechanics techniques can be employed to model the individual constituents within a composite

material. Typically, a repeating unit cell (RUC) in the composite microstructure is identified, and analysis is

performed on that RUC assuming periodic boundary conditions. The response of a point in a continuum is

determined assuming an infinite array of the RUCs. However, representative volume element (RVE) methodologies exist which incorporate applying non-periodic boundary conditions to a subvolume that accurately

represents the composite microstructure. The RVE is meant to represent the actual microstructure of the

continuum, and the size of the features of the microstructure is preserved. Micromechanics can be utilized

to provide the homogenized composite stiffness, or they can be used to model damage and failure within the

constituents. If utilized for the latter, the global mechanisms can arise through the natural evolution and

interaction of the mechanisms in the constituents of the micromechanics model. Numerous micromechanical frameworks exist that encompass analytical, semi-analytical, and numerical techniques. An expansive

review of many micromechanics theories is given in Ref. 1. In this work, a continuum damage model

(CDM) is implemented within the generalized method of cells (GMC), high-fidelity generalized method of

NASA/TM—2012-217649 1

Progressive Failure of a Unidirectional Fiber-Reinforced 

Composite Using the Method of Cells: Discretization 

Objective Computational Results 

Evan J. Pineda and Brett A. Bednarcyk 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Anthony M. Waas 

University of Michigan 

Ann Arbor, Michigan 48109 

Steven M. Arnold 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Abstract 
cells (HFGMC), and finite element method (FEM) to capture the nonlinear response, due to progressive failure, of a unidirectional fiber-reinforced polymer matrix composite (PMC) under a combination of transverse

tension/compression and transverse shear loading. The focus is on developing a semi-analytical methodology

that is insensitive to changes in the density of the numerical discretization, which is necessary for predicting

damage state and failure.

With the method of cells (MOC), developed by Ref. 2, a rectangular composite RUC was discretized

into four subvolumes, called subcells. One of the subcells was occupied by the fiber material and the rest

were occupied by the matrix. Linear displacement fields were assumed in each of the subcells. Displacement

and traction continuity conditions were enforced, in an average integral sense, at the subcell interfaces, along

with periodic boundary conditions at the RUC boundaries to derive a set of equations that yielded a strain

concentration matrix which could, in turn, be used to obtain the local subcell strains from the applied,

global fields. Following determination of the subcell strains, the subcell stresses are readily calculated using

the local constitutive laws, and volume averaging can be used to obtain the homogenized thermomechanical

properties of the composite. MOC was later extended to the generalized method of cells by Ref. 3 which

accommodated any number of subcells and constituents in two-periodic directions. Ref. 4 adapted the

formulation to accommodate triply-periodic materials. Finally, Ref. 5 developed the high-fidelity generalized

method of cells which utilized second order displacement field approximations in the subcells, rather than

linear. Ref. 6 showed the local elastic field accuracy produced by HFGMC corresponded very well to FEM;

whereas, Ref. 7 compared the accuracy of the fields in inelastic phases to composite cylinder assemblage

(CCA). Ref. 8 utilized HFGMC to model fiber-matrix debonding in metal matrix composites, and Ref. 9

implemented a multi-axial damage model in HFGMC. Reformulations, of GMC and HFGMC, which reduced

the total number of unknowns in the problem were introduced by Refs. 10 and 1, respectively. Herein, the

reformulated, doubly-periodic versions of GMC and HFGMC are employed, and the reader is referred to

Refs. 10 and 1 for the complete details on these formulations. These semi-analytical methods offer efficient

alternatives, compared to fully numerical methods such as FEM, for computing local fields within a composite

RUC and are ideal for implementation in a multiscale framework.

The generality of the GMC and HFGMC formulations admit any constitutive behavior at the subcell

level. However if the response of the subcell material exhibits post-peak strain softening, the tangent stiffness

tensor of the subcell loses positive definiteness. This leads to pathologically mesh dependent behavior. This

mesh dependency, if not eliminated, will not provide predictive capability, rather, the capability can be only

used with confidence for “simulations”.

The smeared crack band approach,11 which introduced a characteristic element length into the post-peak

strain softening damage evolution formulation, can be used to alleviate the dependency of failure, and strain

localization, on the discretization size within the continuum domain. The tangent slope of the softening

stress-strain curve was scaled by the characteristic length to ensure that total strain energy release rate

(SERR) upon complete failure (i.e. zero stress) is always equal to the prescribed fracture toughness (or

critical SERR), regardless of the size of the discretization. In the original formulation, the band was always

oriented perpendicular to the direction of maximum principal stress; thus, the crack band always advanced

under pure mode I. Refs. 12 and 13 later reformulated the model to incorporate a fixed crack band that

evolved under mixed-mode conditions. Both formulations employ triangular degradation schemes. Later,

Ref. 14 incorporated more sophisticated initiation criteria to predict the onset of mixed-mode crack bands.

All of these smeared crack formulations assume linear elastic behavior up to the initiation of the crack band,

followed by immediate post-peak strain softening. However, Ref. 15 coupled pre-peak plasticity with crack

band governed post-peak strain softening to model failure of concrete. Recently, a thermodynamically-based

work potential theory was developed which allows for pre-peak progressive damage, as well as, post-peak

progressive failure, governed by the smeared crack band approach, in homogenized laminae,,16 for fiber

reinforced laminates.

A variation of the crack band theory for concrete structures11 is implemented here within GMC, HFGMC

and FEM to model mesh objective failure of continuous fiber-reinforced PMCs. With this implementation,

presented in Section II, two scenarios are considered to determine the mode in which the cracks contained

in the crack band grow. If the principal stress that has the highest magnitude is tensile, it is assumed that

it is more energetically favorable for the crack band to form perpendicular to the maximum principal stress

and for the cracks within the band to advance under mode I conditions. Conversely, if the magnitude of

a compressive principal stress is higher than the other principal stresses, the cracks within the crack band

evolve under mode II conditions (due to internal, Mohr-Coulomb friction) and are oriented with the plane

NASA/TM—2012-217649 2

of maximum shear stress.
Previously, similar studies utilized FEM to model the composite microstructure. Plasticity was employed

to model the non-linear response of the matrix and failure was introduced by allowing the fiber-matrix

interface to debond using cohesive zone elements.17–19 The response of a statistical sample of RUCs subjected

to a combination of transverse compressive and transverse shear was reported. In the current work, the

focus is restricted to the microscale to evaluate the capabilities of the smeared crack band model to predict

progressive failure evolution within a composite microstructure, using the semi-analytical methods (GMC and

HFGMC), by verifying and validating this implementation, as well as the overall utility of the micromechanics

models employed, against available experimental data, and an analogous, fully numerical (FEM) model. The

objectivity of failure evolution with respect to the level of discretization used in HFGMC was previously

demonstrated in Ref. 20. In Section III, the results for an RUC containing 13, randomly placed fibers

subjected to a combination of transverse tension, transverse compression, and transverse shear are presented.

The stress-strain response obtained from GMC and HFGMC are compared to experimental data (where

available) and FEM results. Furthermore, failure path predictions obtained from the three methods are

exhibited.

Using GMC to model an RUC containing 13 randomly arranged fibers is beyond the intended use of this

ultra-efficient, semi-analytical method, and it was expected that with GMC, as with any analytically-based

model, that the results obtained from the simulations of this architecture would not be accurate. These

GMC simulations were executed to give the authors a measure of how the maximum disorder affected the

GMC results. Additional simulations were carried out, utilizing an ordered, square packing architecture

to demonstrate that the accuracy lost with the complex could be regained, and GMC is a viable tool for

modeling progressive failure in the constituents of a composite. Due to the ultra-efficient nature of GMC,

it is far more tractable within a multiscale setting than HFGMC or FEM, and offers fidelity beyond that

available with classical mean field theories. Single fiber RUCs were loaded using the same loading scenarios

used for the multiple fiber RUCs and the results obtained from GMC, HFGMC, and FEM are compared.

The FEM model used to represent the RUC contained a much finer mesh than GMC/HFGMC. This was

to ensure that the FEM solutions were as accurate as possible and provided a fully numerical baseline to

evaluate the performance of GMC and HFGMC. These results are presented in Section IV.

II. Modeling Constituent-Level Post-Peak Strain Softening with the Smeared

Crack Band Approach

GMC and HFGMC are efficient (relative to fully-numerical methods), semi-analytical tools useful for

modeling details of the microstructure of a composite material. Additionally, they are readily amenable for

implementation into a multiscale framework. Although, physics-based, discretization objective, constitutive

theories for modeling progressive failure must be in place to accurately predict the response of a structure

that is failing. For pre-peak loading (i.e. positive-definite tangent stiffness tensor), there are a multitude of

non-linear elastic, plastic, continuum damage mechanics, and visoelastic/plastic theories available that can

predict the evolution of the appropriate mechanisms in the composite. However when the local fields enter

the post-peak regime of the stress-strain laws, most of these theories breakdown in a numerical setting and

display pathological mesh dependence.21, 22

The lack of positive definiteness of the elastic, or inelastic, tangent stiffness tensor leads to imaginary

wave speeds in the material. The longitudinal wave speed in an isotropic material is given by

cL =

s

E(1 − ν)

ρ(1 + ν)(1 − 2ν)

(1)

where cL is the wave velocity, E is the tangent stiffness of the material, ν is the Poisson’s ratio, and ρ is

the material density. A one-dimensional approximation yields v =

qE

ρ

. The existence of an imaginary

wave speed results in a boundary value problem that is ill-posed.21–24 Physically, a material must posses a

positive-definite tangent stiffness tensor, and in fact, at the micro-scale the material tangent stiffness tensor

always remains positive-definite. However for practical purposes, engineers must model structures at scales

much larger than the flaws in the material, and the homogenized continuum representation of a material

containing the nucleation and propagation of discontinuities, such as cracks or voids, exhibits post-peak strain

softening in the macroscopic, homogenized, stress-strain response. This homogenized response is assumed

NASA/TM—2012-217649 3
to govern over a suitable volume of the material, appropriate to the microstructure of the material. Loss

of positive-definiteness of the tangent stiffness tensor leads to a material instability, which manifests as a

localization of damage into the smallest length scale in the continuum problem.24 In GMC, or HFGMC, this

is a single subcell. Thus, the post-peak softening strain energy is dissipated over the volume of the subcell

that the damage localizes to. Since a stress-strain relationship prescribes the energy density dissipated during

the failure process, the total amount of energy dissipated in the subcell decreases as the size of the subcells

is reduced, and in the limit, zero energy is required to fail the structure. An illustration showing this size

dependence is given in Figure 1.

A simple way to remedy this non-physical behavior in a numerical setting is to judiciously scale the

post-peak softening slope of the stress-strain constitutive law. Then, the failure energy density dissipated

becomes a function of the characteristic length of discretized continuum. Refs. 25 and 11 first proposed

a crack band model in which post-peak softening damage (herein referred to as failure) in the material

was assumed to occur within a band. The post-peak slope of the material constitutive law was scaled by

the characteristic length of the finite element exhibiting failure; such that, the total SERR in the element,

upon reaching a state of zero stress, and the material fracture toughness were coincident. In this reference,

equivalence between this smeared crack approach and a line crack approach is presented. Subsequently,

Ref. 26 exhibited propagation of a crack band not aligned with the mesh bias. In this work, the crack

band theory is implemented within the GMC and HFGMC micromechanics models, in the MAC/GMC suite

of micromechanics codes developed at the NASA Glenn Research Center,27, 28 and used to analyze crack

band growth in composite RUCs. The results of the semi-analytical methods are compared to experimental

data, where available, and an equivalent FEM models containing the same crack band implementation. The

following subsections provide theoretical details on the crack band model.

II.A. Physical Behavior of Crack band

The smeared crack band model is meant to capture the behavior of a region of a material wherein numerous

microcracks have initiated, and they coalesce to form a larger crack. Figure 2 displays a crack band of width

wc embedded in a continuum. The domain of the crack band is denoted as Ω0 and the remaining continuum

as Ω. The crack band is oriented within the continuum such that, for a given point within the crack band,

the unit vector normal to the crack band is n.

The total energy dissipated during the failure process is dissipated over Ω0, and the size wc of Ω0is a

material property directly related to the material fracture toughness.11

wc =

2GC

σ

2

C

μ

1

E

−

1

ET

¶−1

(2)

where σC is the critical stress for initiation of the post-peak regime in the 1D material stress-strain law

(given in Figure 3a), and ET is the negative tangent slope in that regime. The fracture toughness GC , or

critical strain energy release rate, of the material is given by the area under the 1D traction-separation law

(Figure 3b) which governs the cohesive response of the separation of crack faces as a crack propagates in the

material. The energy density dissipated during failure WF is related to the material fracture toughness by

the characteristic length in the material.

GC = wcWF (3)

II.B. Implementation of Crack Band Model in GMC/HFGMC

The crack band model is implemented in the GMC and HFGMC micromechanical framework. The local

subcell fields are used to govern crack band evolution in the constituents of the composite. Figure 4 shows

the discretization of the continuum displayed in Figure 2. A magnified view of the crack band embedded in

a single subcell is also displayed in Figure 4. Since all of the energy dissipated in the crack band is smeared

over the subcell volume, the subcell must be large enough to contain the crack band of width wC . Note that

Figure 4 shows a 2D geometry for illustrative purposes, but the crack bands can also evolve in a general 3D

space.

NASA/TM—2012-217649 4
II.B.1. Mode I Crack Band Growth Under Tensile Principal Stress

The orientation of the crack band in subcell βγ is given by the vector n

(βγ)

1

(see Figure 4) and is determined

from the local principal stress state ( ̄σ

(βγ)

1

, σ ̄

(βγ)

2

, σ ̄

(βγ)

3

). In a monolithic material, cracks orient such that

the crack tips are always subjected to pure mode I (opening mode) conditions unless there are constraints

that limit the crack orientation. In a micromechanical analysis the composite material is composed of

separate, monolithic constituents; thus, in the matrix, the crack band runs perpendicular to  ̄σ

(βγ)

1

, the

principal stress with the largest magnitude, |σ ̄

(βγ)

1

| > |σ ̄

(βγ)

2

| > |σ ̄

(βγ)

3

|, if  ̄σ

(βγ)

1 ≥ 0 (tensile). Under these

conditions, a crack oriented as such, is subjected to pure mode I loading, locally. Although, the resulting

global behavior may appear to be mixed mode because of the influence of the fibers on the matrix crack

band path. Crack orientation and evolution is determined differently if  ̄σ

(βγ)

1 < 0 as described in Section

II.B.2. The characteristic length of the subcell l

(βγ)

C

is determined as the dimension of the subcell running

parallel to n

(βγ)

1

. Crack band initiation is determined using a very simple, but physical,29 maximum stress

criterion.

σ ̄

(βγ)

1

σ

(βγ)

C

= 1, σ ̄

(βγ)

1 ≥ 0 (4)

where σ

(βγ)

C

is the cohesive strength of the crack band. Once the crack band has initiated, the crack band

orientation is fixed as time evolves.

Once the crack band orientation has been calculated, the subcell compliance is rotated into the principal

frame using the transformation matrix.

T = [n

(βγ)

1 n

(βγ)

2 n

(βγ)

3

][e1e2e3] (5)

where n

(βγ)

1

, n

(βγ)

2

, and n

(βγ)

3

are the principal stress directions, and e1, e2, and e3 are the unit basis

vectors. All material degradation due to crack band evolution is imposed on the rotated compliance S ̄(βγ),

the components of which are given by:

S ̄

(βγ)

ijkl = TpiTqjS

(βγ)

pqrsTkrTls (6)

The strain energy released during the formation of new surfaces corresponding to the growth of cracks

within the crack band is assumed to be dissipated over the entire subcell volume. Therefore, the post-peak

softening slope E

(βγ)

IT , and the strain at which a the principal stress state is zero, is calculated using the

characteristic length of the subcell l

(βγ)

C

and the material fracture toughness G

(βγ)

IC (see Figure 3a).

2

(βγ)

F =

2G

(βγ)

IC

σ

(βγ)

C

l

(βγ)

C

(7)

E

(βγ)

IT =

Ã

1

E ̄

110

−

2

(βγ)

F

σ

(βγ)

C

!−1

(8)

where E ̄

110 is the undamaged, axial Young’s modulus in the principal frame. It should be noted that E

(βγ)

IT

must be less than zero; therefore, by Eqs. (7) and (8), a restriction is placed on the maximum allowable

subcell size.

l

(βγ)

C <

2G

(βγ)

IC E ̄

110

σ

(βγ)

C

(9)

The local, rotated, subcell strain state  ̄2

(βγ)

i







2 ̄

(βγ)

1

2 ̄

(βγ)

2

2 ̄

(βγ)

3







=







T

(βγ)

1i

2 ̄

(βγ)

ij T

(βγ)

1j

T

(βγ)

2i

2 ̄

(βγ)

ij T

(βγ)

2j

T

(βγ)

3i

2 ̄

(βγ)

ij T

(βγ)

3j







(10)

NASA/TM—2012-217649 5
is used to degrade the rotated compliance components. The scalar, mode I, damage factor D

(βγ)

I

is calculated

using the rotated strain corresponding to  ̄σ

(βγ)

1

.

D

(βγ)

I = 1 +

E

(βγ)

IT 3

2

(βγ)

C − 2 ̄

(βγ)

1

 ́

E ̄

1102 ̄

(βγ)

1

(11)

where 2

(βγ)

C

is the value of  ̄2

(βγ)

1 when the initiation criterion, Eq. (4), is satisfied. If D

(βγ)

I

is less than zero,

no damage occurs, and a maximum damage level of one corresponds to a zero stress state on the softening

stress-strain curve. Also, damage healing is inadmissable.

D ̇

(βγ)

I ≥ 0 (12)

Components of the rotated compliance matrix are degraded with the damage factor.

S ̄(βγ) =































































S ̄0

(βγ)

1111

3

1 − D

(βγ)

I

 ́ S ̄

(βγ)

1122 S ̄

(βγ)

1133 0 0 0

S ̄

(βγ)

1122 S ̄

(βγ)

2222 S ̄

(βγ)

2233 0 0 0

S ̄

(βγ)

1133 S ̄

(βγ)

2233 S ̄

(βγ)

3333 0 0 0

0 0 0 S ̄

(βγ)

2323 0 0

0 0 0 0

S ̄0

(βγ)

1313

3

1 − D

(βγ)

I

 ́ 0

0 0 0 0 0

S ̄0

(βγ)

1212

3

1 − D

(βγ)

I

 ́































































(13)

Since the crack band orientation is fixed upon initiation, the S ̄

(βγ)

1313 and S ̄

(βγ)

1212 shear compliances in the rotated

frame are degraded, as well as the S ̄

(βγ)

1111 compliance, so that the crack band faces normal to n

(βγ)

1

are free of

normal and shear tractions when all of the crack band energy has been dissipated (i.e. l

(βγ)

C W

(βγ)

F = G

(βγ)

IC ).

Once the compliance in the rotated frame is degraded, the compliance is transformed back to the global

frame to yield the new subcell compliance.

S

(βγ)

ijkl = T

−1

pi T

−1

qj S ̄(βγ)

pqrsT

−1

kr T

−1

ls (14)

Note that, damage introduced in the principal frame, through Eq. (13), can induce normal-shear coupling

in the global frame.

II.B.2. Mode II Crack Band Growth Under Compressive Principal Stress

Crack band growth under pure mode I conditions are energetically favorable; however, conditions arise, under

which, mode I crack growth is not possible. For instance, the crack cannot grow under in mode I if the normal

traction at the crack tip is compressive. This occurs when the principal stress with the maximum magnitude

is compressive ( ̄σ

(βγ)

1 < 0). However, experimental data shows that monolithic materials subjected to global

compression will eventually fracture.

In brittle and quasi-brittle materials, it has been shown that local, internal friction results in shear (mode

II) fracture when the monolithic material is subjected to compressive stresses.30–33 Consequently, if  ̄σ

(βγ)

1 < 0

it is assumed that the crack band is aligned with the plane of maximum shear stress  ̄τ

(βγ)

in the matrix

NASA/TM—2012-217649 6
subcell, and a Mohr-Coulomb failure criterion will be used to indicate crack band initiation under maximum

principal compression.

τ ̄

(βγ)

E

τ

(βγ)

C

= 1, σ ̄

(βγ)

1 < 0 (15)

where τ

(βγ)

C

is the cohesive shear strength of the matrix, and  ̄τ

(βγ)

E is an effective shear stress that includes

the influence of the traction normal to the crack band  ̄σ

(βγ)

n .

τ ̄

(βγ)

E = |τ ̄

(βγ)

| + μiσ ̄

(βγ)

n

(16)

where μiis the internal friction coefficient and must lie between 0 < μi ≤ 1.5 to obtain physically reasonable

surface friction coefficient.33 Since the principal stress with the largest magnitude is compressive, the traction

normal to the maximum shear stress plane  ̄σ

(βγ)

n must also be compressive. Hence, in Eq. (16), an increase

in the normal compressive stress will yield a reduction in the effective shear stress and an increase in

apparent shear strength.34 derived relationships between internal friction coefficient, shear strength, the

tensile strength and compressive strength of the material, yielding an expression for μiin terms of σ

(βγ)

C

and

τ

(βγ)

C

.

μi = tan







sin−1







3

2τ

(βγ)

C

 ́2

− σ

(βγ)

C

2

¡

2 ̄τ

(βγ)

¢2

+ σ

(βγ)

C

2













(17)

The Mohr-Coulomb theory has been used as yield criterion for the matrix of a PMC, experiencing local

compression, in an FEM micromechanics model.17, 18 Additionally, similar criteria are used for shear failure

under compressive loading in homogenized composite materials.35–37 In these theories, the orientation of

the crack is not aligned with the plane of maximum shear stress because of the influence of the fibers on the

crack path, but rather, orientation is an input to the theories obtained from experimental data. Since the

influence of the fibers are explicitly accounted for in the present methodology, the crack may remain oriented

with the maximum shear stress plane, and the local stress in the matrix will drive the crack path.

Subsequent to mode II crack band initiation via Eq. (15), the compliance tensor of the subcell is rotated

into the maximum shear stress frame using Eqs. (5) and (6), where n

(βγ)

1

is a unit-vector perpendicular to

the plane of maximum shear stress. Features of the post-peak softening shear stress-shear strain curve, such

as the endpoint and post-peak tangent stiffness, can be calculated from the mode II fracture parameters.

γ

(βγ)

F =

2G

(βγ)

IC

τ

(βγ)

C

l

(βγ)

C

(18)

E

(βγ)

IIT =

Ã

1

G ̄

120

−

γ

(βγ)

F

τ

(βγ)

C

!−1

(19)

where G ̄

120 is the undamaged, axial shear modulus in the maximum shear stress frame. Again the tangent

stiffness E

(βγ)

IIT must be less than zero; therefore, Eqs. (18) and (19) place a restriction is on the maximum

allowable subcell size.

l

(βγ)

C <

2G

(βγ)

IC G ̄

120

τ

(βγ)

C

(20)

The shear strain corresponding to the maximum shear stress  ̄γ

(βγ)

is obtained by rotating the strain

tensor.

γ ̄

(βγ) = T

(βγ)

1i

2 ̄

(βγ)

ij T

(βγ)

2j

(21)

An effective shear strain  ̄γ

(βγ)

E

can be defined that is work conjugate with the effective shear stress  ̄τ

(βγ)

E

.

τ ̄

(βγ)

E

dγ ̄

(βγ)

E =  ̄τ

(βγ)

dγ ̄

(βγ) +  ̄σ(βγ)

n d2 ̄

(βγ)

n (22)

where  ̄τ

(βγ) and  ̄σ

(βγ)

n are the shear and normal tractions acting on the crack faces oriented parallel to the

plane of maximum shear stress, and  ̄γ

(βγ) and  ̄2

(βγ)

n are the apparent shear and normal strains with respect

NASA/TM—2012-217649 7
to the maximum shear orientation in the subcell, including the effects of crack tip opening displacement.

Taking the derivative of Eq. (22) with respect to  ̄τ

(βγ) yields

∂τ ̄

(βγ)

E

∂τ ̄

(βγ)

dγ ̄

(βγ)

E = dγ ̄

(βγ)

(23)

and along with the derivative of Eq. (16) with respect to  ̄τ

(βγ)

∂τ ̄

(βγ)

E

∂τ ̄

(βγ)

= 1 (24)

finally, after integrating, gives:

γ ̄

(βγ)

E =  ̄γ

(βγ)

(25)

which states that the shear strain in the rotated, maximum shear stress frame and the effective shear strain

that is work conjugate to the effective, Mohr-Coulomb shear stress are equivalent.

Thus,  ̄γ

(βγ)

is used to degrade the rotated compliance components. The mode II damage factor D

(βγ)

II is

calculated using  ̄γ.

D

(βγ)

II = 1 +

E

(βγ)

IIT 3

γ

(βγ)

C − γ ̄

(βγ)

 ́

G ̄

120γ ̄

(βγ)

(26)

where γ

(βγ)

C

is the value of  ̄γ

(βγ) when the initiation criterion, Eq. (15), is satisfied. In the case of mode II

fracture, only the shear moduli are degraded. The normal direction is subjected to compression, and thus,

retains its stiffness.

S ̄(βγ) =

























































S ̄

(βγ)

1111 S ̄

(βγ)

1122 S ̄

(βγ)

1133 0 0 0

S ̄

(βγ)

1122 S ̄

(βγ)

2222 S ̄

(βγ)

2233 0 0 0

S ̄

(βγ)

1133 S ̄

(βγ)

2233 S ̄

(βγ)

3333 0 0 0

0 0 0 S ̄

(βγ)

2323 0 0

0 0 0 0

S ̄0

(βγ)

1313

3

1 − D

(βγ)

II  ́ 0

0 0 0 0 0

S ̄0

(βγ)

1212

3

1 − D

(βγ)

II  ́

























































(27)

Then, the compliance tensor is rotated back to the global frame using Eq. (14). It is possible, due to

load redistribution within the composite microstructure, that the normal component (in the maximum shear

frame) is tensile following mode II crack band initiation via the Mohr-Coulomb criterion. Accordingly, the

S ̄

(βγ)

1111 component of the compliance matrix is abruptly degraded using the previously attained damage state

if this situation arises.

S ̄

(βγ)

1111 =







S ̄0

(βγ)

1111

3

1 − D

(βγ)

II  ́ σ ̄

(βγ)

n > 0

S ̄0

(βγ)

1111 σ ̄

(βγ)

n ≤ 0

(28)

III. Modeling Crack Band Growth in RUC Containing 13 Randomly

Arranged Fibers

III.A. Details of Numerical Models

The smeared crack band theory, presented in Section II.B, was implemented within GMC, HFGMC, and

FEM, and it was used to model the progression of post-peak softening failure within a unidirectional fiberNASA/TM—2012-217649 8
reinforced composite RUC. The main objective was to capture the failure evolution of a composite RUC

under transverse tension, transverse compression, transverse shear, and a combination of transverse normal

and shear loading. In tension, the lamina is susceptible to transverse cracking, as evidenced by the scanning

electron microscope (SEM) image of a unidirectional composite failed under tension in Figure 5. The behavior

of this mechanism is observed to be quite brittle, as is evidenced by tension tests on 90◦coupon laminates.

However, in compression a network of shear bands develop, accompanied by some matrix cracking and fibermatrix debonding (see the SEM image in Figure 6), and the stress-strain response is much more ductile.

To emulate these failure modes, a doubly-periodic RUC of the 2-3 plane of the composite was created (see

Figure 7). This RUC was subjected to global normal and shear strains. The constitutive behavior of the

matrix subcells follows linear elasticity, coupled with the crack band model for post-peak strain softening

formulated earlier. This model is verified against experimental data.

A representation of the 2-3 plane of a typical composite lamina (where the 1-axis runs in the longitudinal

direction of the fibers, and the 2- and 3-axes span the plane of transverse isotropy) was produced to examine

the evolution of transverse cracks within the composite. The blueprint for the discrete model was created

by randomly placing thirteen circular fibers with a diameter of 5 μm in a square box. The domain was

then discretized into a number of subcells. The considered formulations of GMC and HFGMC only admit

rectangular subcells in the doubly-periodic models (parallelepiped in triply-periodic models); so, some of

the fiber domains would overlap or lie directly adjacent to one another (an isoparametric formulation of

HFGMC has been recently developed by Ref. 38). Mesh objective, subcell interfacial debonding has been

formulated for GMC (see Ref. 39, 40) but has not been fully developed within MAC/GMC. As a result, thin

matrix channels were inserted between any adjacent fibers to avoid arrest of the crack band for non-physical

reasons such as the inability of two adjacent fiber to separate. The final architecture and discretization is

shown in Figure 7, which contains 81 subcells x 85 subcells. The blue subcells indicate fiber material and

the green represent the matrix. The dimensions of the RUC (21.25 μm x 21.25 μm) were chosen such that

the final fiber volume fraction in Figure 7 would be preserved at 58%, corresponding to the experimental

data provided in Ref. 41. The RUCs were loaded in strain control,and MAC/GMC was used to perform the

calculations.

The original intention of this study was to verify and validate the use of the crack band theory within a

micromechanics model. Thus, a complex RUC (which is intractable in a multiscale setting) was chosen to

demonstrate its utility. This highly disordered architecture is beyond the recognized capabilities of GMC,42

and it was expected that the GMC solutions would contain significant error. However, the GMC simulations

of this architecture were executed to discern the performance of GMC given an extreme case. The following

section forgoes the complicated RUC in lieu of an ordered structure, for which, GMC is more applicable.

The elastic properties for a Silenka E-glass/MY750/HY917/DY063 lamina were taken from data from the

world wide failure exercise WWFE.41 Elastic properties for Silenka E-glass fiber and MY750/HY917/DY063

epoxy matrix used are given in Table 1. Note that the matrix properties were correlated to represent the

in-situ properties of the matrix, which differ significantly from those of the neat material,43 and match the

global properties of the composite. As a result of this calibration, the isotropic relationship between the

shear modulus and the Young’s modulus and Poisson’s ratio of the matrix is not maintained exactly. The

fracture properties (also Table 1) were calibrated to produce and ultimate transverse tensile stress of 40

MPa, and the compressive properties were calibrated to match transverse compression data, reported in Ref.

41.

The same properties and mesh were also used in an analogous FEM model. The domain was discretized

using 2D, quadrilateral, generalized plane strain, reduced integration CPEG4R elements in the commercially

available Abaqus finite element software.44 Periodic boundary conditions were imposed on the RUC using

the *EQUATION keyword, which ties the displacement of corresponding node pairs to a global displacement.

u

1

i − u

2

i =  ̄u

R

i

(29)

Nodes on the left and right edges were tied to a reference node using Eq. (29), while nodes on the top and

bottom edges were tied to an additional reference node. Displacements were applied to the reference nodes

to match the global applied strains in the GMC/HFGMC models. This procedure is similar to that used in

Ref. 18 to apply periodic boundary conditions. Additionally, the out-of-plane displacement of all nodes was

tied to a reference node to enforce the generalized plane strain condition. This node was allowed to move

freely in 1-direction while the out of plane rotations were constrained ensuring that the front and back plane

of the RUC remained parallel throughout the simulation. An implicit, dynamic solver was used to solve the

NASA/TM—2012-217649 9
equations of motion as the RUCs were loaded, to the same global strain levels as the GMC and HFGMC

models, over a total time 1000s. An apparent density of 1.57 g/cm3 was assigned to the fiber and matrix

constituents to yield reasonable masses. The smeared crack band constitutive model in Section II.B was

implemented in FEM using an Abaqus/Standard UMAT User Subroutine44 for the matrix elements. The

fiber was treated as linear elastic.

III.B. Results: Transverse Tension

Figure 8 shows the transverse stress (σ22) versus transverse strain (222) response obtained from applying

a global transverse strain to an RUC containing 13 randomly arranged fiber (see Figure 7) using GMC,

HFGMC and FEM to model the RUC. The tensile fracture properties (σ

(βγ)

C

and G

(βγ)

IC ) used in all three

models were obtained by calibrating them such that the ultimate, global transverse stress from the HFGMC

model correlated with the documented ultimate transverse stress from Ref. 41, which was 40 MPa. Only a

single data point given in Ref. 41, and it was reported that the tensile behavior was linear up to failure.

The correlated, ultimate, global transverse stress computed using HFGMC was 40.2 MPa at a strain

of 0.00248, and the material response was linear up to failure. After the ultimate stress was achieved, the

global stress dropped suddenly. However, the RUC retained some load carrying capability, and the global

stress increased with applied global strain for a short period until the global stress reached 34.1 MPa at an

applied strain of 0.00264. After which, the load carrying capability of the RUC was completely diminished.

When the same fracture properties were used in an analogous FEM model, the predicted ultimate stress was

41.7 MPa at 222 = 0.00253, 3.7% higher than the HFGMC results. The FEM results were also linear up

to ultimate failure. Although, the FEM results exhibited two additional stress peaks at 35.8 MPa and 36.3

MPa corresponding to applied strain levels of 0.0028 and 0.00324, respectively. Both HFGMC and FEM

provided comparable quantitative results for applied, tensile transverse strain.

As expected, the results obtained from GMC differ significantly from both HFGMC and FEM. The initial

slope of the σ22-222 curve produced by GMC was lower than the slope calculated by HFGMC and FEM.

The first drop in the global stress occurred at σ22 = 41.5 MPa (222 = 0.0032); however, the ultimate stress

was considerably higher (σ22 = 48.2 MPa, 222 = 0.004). Additionally, the post-peak response obtained from

GMC exhibited many more stress peaks, and less severe stress drops than HFGMC and GMC. Furthermore,

the overall transverse response of the RUC was very tough using GMC, whereas, HFGMC and FEM showed

very brittle behavior. The discrepancy between GMC and HFGMC (or FEM) is attributed to the lack of

shear coupling, inherent in the GMC (and any other first order theory) formulation. Closely packed fibers

would promote an increase in the shear stresses between those fibers while the shear stress in matrix rich

region would be minimal. Since the global loading is strictly normal, and there is no shear coupling in GMC,

the shear stress calculated by GMC is zero everywhere. This yields an error in the overall equilibrium of

the structure and results in a lower stiffness and erroneous estimate of the ultimate strength of the RUC.

However, the random arrangement of the fibers in the RUC amplifies this drawback of GMC. This detail in

this architecture exceeds the applicable scope of GMC, as GMC promotes solution efficiency over complexity.

Furthermore, the performance of GMC is vastly improved when the microstructure of the RUC is ordered,42

as is demonstrated in the following section.

Contours of the mode I degradation variable are obtained from transverse tensile loading of an RUC using

HFGMC, FEM, and GMC are presented in Figures 9-11, respectively. Using HFGMC (see Figure 9), a small

crack band developed between two fibers, and can be seen in Figure 9a (note that fibers are not displayed

in the HFGMC contour plots), when the RUC reached the ultimate stress of 40.2 MPa. At the second peak

stress, given in Figure 9b, (σ22 = 34.5 MPa), the original crack band progressed in the negative x3-direction

and in the positive x3 direction across the periodic boundary. Also, an additional crack band began to form

elsewhere in the RUC. Finally, once the load carrying capability of the RUC completely dropped, the tips of

the first crack band navigated around the fiber inclusions and met, forming a fully bridged crack band that

spans across the x3-dimension of the RUC, as shown in Figure 9c.

A similar crack band evolution was apparent in the FEM results exhibited in Figure 10. The development

of the crack band at the ultimate stress (Figure 10a), and the initial progression of the crack band (Figure

10b), was similar to what was observed with HFGMC. However, the discrepancy in the σ22-222 response

between HFGMC and FEM, subsequent to the second peak stress, can be attributed to the contrast in

the evolution of the crack band between the two models. With the FEM model, the RUC reloaded after

the second stress peak was reached. This corresponded to the stagnation of the initial crack band, but the

growth of a second crack band (see Figure 10c). In the HFGMC model, the first crack band continued to

NASA/TM—2012-217649 10
grow and resulted in a complete reduction of the global transverse stress after the second stress peak was

achieved. Finally, in the FEM model, a third stress peak was observed (Figure 8), after which, the load

carrying capability of the FEM RUC was severely reduced and the first crack band linked up with the second

crack band forming a single crack band that progressed along the x3-direction of the RUC. Note that the

crack band path predicted by HFGMC and FEM match well with the crack paths observed in experiments,

shown in Figures 5

Figure 11 displays the mode I crack band evolution calculated using GMC as the micromechanics model,

and as presumed, exhibited very different behavior than HFGMC or FEM. Since the uniform subcell strains

are coupled through displacement continuity conditions governing GMC,3strain, and hence damage, localization cannot occur within the a complex RUC. The initiation of the crack band was similar to HFGMC

and FEM (see Figure 11a), however; the progression of crack bands subsequent to the ultimate stress was

tremendously sporadic (see Figures 11b-11d), and the a fully bridged crack band never developed. The large

amount of energy released during the transverse failure of the GMC RUC can be attributed to this speckled

progression of crack bands throughout the matrix in the RUC.

III.C. Results: Transverse Compression

The transverse stress versus transverse strain results for RUCs containing 13 randomly arranged fibers

subjected applied transverse compression are presented in Figure 12. For this loading case, the mode II

strength (τ

(βγ)

C

) and critical SERR (G

(βγ)

IC ) were calibrated within HFGMC to match the experimental data41

in Figure 12 as closely as possible. The transverse compressive stress in the HFGMC model dropped, and

the RUC failed, at σ22 = -99.8 MPa, far before the reaching the end of the experimental data. The end

of the FEM data at σ22 = -108.4 MPa corresponded with divergence of the FEM solution, indicating that

the solution was near failure as well. As with the tensile results, the transverse stiffness and strength (-62.4

MPa) of the GMC model under compression did not correspond to the HFGMC or FEM models.

Figure 13 exhibits the failure contours obtained from the HFGMC simulations. Figures 13a-13c show

the mode II failure contours, resulting from local principal compression. Figures 13d-13e show the mode I

failure contours, resulting from local principal tension. Finally, the superposition of the mode I and mode

II failure contours are displayed in Figures 13g-13i. Prior to the ultimate stress, an array of mode II crack

bands developed throughout the RUC (Figure 13a). These mode II crack bands were accompanied by three

mode I crack bands, shown in Figure 13d. At the ultimate stress there was a large degree of mode II damage

(Figure 13b), and the most prevalent mode I crack band developed fully (Figure 13c). It can be seen in the

superposition plot at the ultimate stress, Figure 13h, that the mode I crack band bridged two separate mode

II crack bands. Upon final failure (Figures 13c, 13f, and 13i), several mode II crack bands were linked and

bridged by a mode I crack band. The mode I crack band was present in the matrix subcells at the boundary

of a fiber in the RUC. This failure mode was consistent with the observed, experimental compressive failure

displayed in Figure 6, which shows numerous shear bands being bridged by tensile fiber-matrix debonding.

The failure contours obtained from the FEM simulations of an RUC under transverse compression are

presented in Figure 14. Prior to the end of the simulation (which corresponded to divergence in the FEM

solution), a small amount of dispersed mode II crack bands was present (see Figure 14a). Additionally, some

mode I crack bands started developing (Figure 14c). At the final time increment, more mode II crack bands

evolved (see Figure 14b), and some mode I crack bands fully developed, as shown in Figure 14d. As with

HFGMC, the mode I crack band observed in the FEM solutions surrounded a fiber in the RUC and bridged

multiple mode II crack bands. This is displayed in Figure 14f, which displays the superposition of mode I

and mode II failure. Unfortunately, the FEM solution diverged before a fully bridged crack could be formed.

The GMC simulation of an RUC subjected to transverse compression exhibited only mode II failure,

which is presented in Figure 15. As with the tensile case, the failure pattern was highly irregular, exhibiting

no localization, and did not represent any physical failure mode.

Since the local mode II fracture toughness (G

(βγ)

IC ) governing local compressive failure was significantly

higher than the mode I fracture toughness (G

(βγ)

IC ), it was expected that the response of the RUC under

transverse compression would be far tougher than under tension. Some toughness was observed, but the

HFGMC and FEM RUCs failed far before the end of the experimental data. Increasing G

(βγ)

IC did not delay

the final failure of the RUC but simply increased the nonlinear slope of the transverse stress-strain curve

after crack band initiation. The presence of fully developed mode I crack band in HFGMC and FEM RUCs

subjected to transverse compression at the ultimate compressive stress led to the suspicion that the mode I

NASA/TM—2012-217649 11
crack bands were driving the final failure of the RUCs, since the final failure was insensitive to the mode II

fracture toughness. Additional simulations were executed were mode I failure was disabled, and degradation

of S ̄

(βγ)

1111, see Eq. (28), was prohibited even if the stress normal to the plane of maximum shear stress was

tensile.

The transverse stress versus transverse strain response of the HFGMC, GMC, and FEM RUCs are plotted

in Figure 16. Both HFGMC and FEM exhibit a much tougher response than displayed in Figure 12, and are

reasonably close to each other. The small discrepancy between the two models indicates a slight difference

in the local maximum shear stresses. Additionally, an ultimate stress is not achieved within a global applied

transverse strain of -0.014. This indicates that the mode II crack band evolution is governing the nonlinearity

in the transverse compressive stress-strain response, and mode I crack band formation is dictating the final

failure. Furthermore, it indicates that thirteen fibers may not be enough represent the actual unit cell within

a composite experiencing strain localization under compression. In a larger RUC, the effects of the formation

of mode I crack band may not be as significant. Further studies will be performed to evaluate the effect of

the size of the RUC on the predicted response of the composite. The response of the GMC RUC is similar to

the simulation that allowed for mode I crack bands and degradation of S ̄

(βγ)

1111 in the maximum shear frame

as expected, since no mode I failure arose in that scenario. Figures 17 and 18 show the mode II crack band

evolution in the HFGMC and FEM RUCs, respectively. Both models exhibited a similar crack band pattern;

however, the crack bands in the HFGMC model are more developed which resulted in the softer stress-strain

behavior.

III.D. Results: Transverse Shear

In the previous two sections, the crack band inputs were calibrated (using HFGMC) to match experimental

data. The same inputs were used in HFGMC, GMC, and FEM, and the performance of the models were

compared. Ref. 41 also provided data for a composite subjected to axial shear (τ12-γ12). However, the

evolution of shear microcracking is heavily influenced by the 3D effects of the composite lamina.43 Adjacent

fibers act to arrest angular microcracks, which leads to a reasonably slow progression of microcracks throughout the lamina. Moreover, constraining plies in a [±45◦

]S laminate, that is typically tested to measure the

axial shear behavior of a composite, add further stability to the damage progression. To model this using

inputs that were calibrated using an unconstrained RUC with transverse tensile and compressive data, the

geometry would need to be accounted for correctly using a fully 3D model. This simulation is outside the

scope of this paper, and is left for future work.

However, more complicated loading scenarios were desired to validate the use of the crack band theory

within the HFGMC micromechanics model. Therefore, transverse shear loading cases were performed and

FEM results are used as the baseline. Figure 19 shows the global transverse shear stress versus transverse

shear strain (τ23-γ23) response of the RUCs subjected to an global applied transverse shear strain. There

was a slight discrepancy between the slopes obtained from HFGMC and FEM. Additionally, the FEM model

exhibited a slight drop in the global transverse shear stress, at γ23 = 0.0043 and τ23 = 28.1 MPa, prior to

attaining the ultimate stress. However, the ultimate transverse shear stresses predicted by the two models

were comparable with HFGMC predicting failure at τ23 = 33.4 MPa, and FEM predicting failure at τ23 =

33.3 MPa. The global shear stiffness calculated using GMC was well below that from HFGMC and FEM,

and the ultimate shear stress was predicted to be 49.2 MPa.

No mode II failure was observed in any of the models. The mode I crack band progressions obtained

from HFGMC and FEM are given in Figure 20 and 21. At the ultimate shear stress a single crack band

forms in both HFGMC (Figure 20b) and FEM (Figure 21b); however, the location of those crack bands

are different. The final failure pattern predicted by HFGMC is shown in Figure 20c, and shows a mode I

crack band progressing at an angle through the entire RUC. A similar failure pattern was predicted by FEM

(Figure 21c); however, the FEM solution diverged before the crack band could fully develop and the global

shear stress could drop to a trivial value. Finally, Figure 22 shows the final mode I degradation pattern

predicted by GMC. A very small level of local failure yielded global failure of the specimen.

III.E. Results: Combined Transverse Normal-Transverse Shear Loading

Figure 23 shows the predicted σ22-222 and τ23-γ23 responses of HFGMC, GMC, and FEM multiple fiber

RUCs subjected to applied transverse tensile strain and transverse shear strain with a ratio of normal to

shear strain of 0.02. Initially the slopes of the normal stress-strain curves predicted by HFGMC and FEM

NASA/TM—2012-217649 12
were the same, whereas, the slope of the shear stress-shear strain curve obtained from FEM was slightly higher

than HFGMC. Also, there was a slight drop in the shear stress predicted by FEM before the ultimate shear

stress was achieved. Both HFGMC and FEM predict similar ultimate transverse tensile and transverse shear

stresses. Additionally, subsequent to the ultimate normal stress, the global normal stress became compressive

in the HFGMC and FEM models. The post peak shear responses of HFGMC and FEM differed, in that

the post peak stress achieved by HFGMC was higher than FEM. Both the stiffnesses and ultimate stress

predicted by GMC were different than HFGMC and FEM. Furthermore, the post-peak shear response of the

GMC RUC was vastly different from HFGMC and FEM.

No mode II failure was present in either the HFGMC or FEM RUCs under the combined transverse

tension-transverse shear loading scenario. The mode I failure patterns for HFGMC and FEM are presented

in Figure 24 and 25. The HFGMC simulation (Figure 24) predicted the initiation of two crack bands: one in

the middle of the RUC and one near the periodic boundary of the RUC. As these crack bands evolve, they

eventually bridged into a single crack that spanned the entire RUC at an angle. The FEM model, shown in

Figure 25, predicted a different location for the crack band initiation and the final crack band (Figure 25c)

progressed in a similar manner to the crack band in the HFGMC model. GMC predicted a minimal amount

of mode I and mode II degradation, as evidenced by Figure 26.

The σ22-222 and τ23-γ23 behavior of an RUC subjected to applied transverse compressive strain and

transverse shear strain, at a ratio of transverse strain to shear strain of -0.75, predicted by HFGMC, GMC,

and FEM is plotted in Figure 27. The FEM solution diverged before a clear ultimate transverse stress

could be obtained. However, prior to this, the FEM and HFGMC models exhibited similar transverse

behavior. The ultimate shear stresses predicted by HFGMC and FEM are comparable, and the post-peak

shear responses are similar. GMC underpredicted the ultimate transverse and shear stresses and exhibited

much softer transverse and shear behavior.

Mode I and mode II degradation progression predicted with HFGMC was displayed in Figure 28 The

response is similar to the loading case incorporating only transverse compressive strain loading; however,

more mode I crack bands were observed that bridge the network of mode II crack bands throughout the RUC.

Analogous failure contours from the FEM model, prior to divergence of the FEM solution, are presented

in Figure 29. As with HFGMC, numerous mode I crack bands are observed bridging mode II crack bands.

However, these crack bands did not develop fully before the FEM solution diverged. GMC exhibited no mode

I failure, but the mode II failure contours are presented in Figure 33. As with the other loading scenarios,

the failure pattern predicted by GMC was disperse and non-physical.

IV. Modeling Crack Band Growth in RUC Representing Square Packing

Architecture

IV.A. Details of Numerical Models

In the previous section, the implementation of the crack band model within HFGMC was verified to capture

the failure response of a unidirectional E-glass/Epoxy composite against experimental data and was shown

to provide results comparable to an FEM model with the same crack band implementation. In this section

the performance of simpler HFGMC and GMC models, more amenable to multiscale implementations, are

compared directly compared to a similar FEM model. Additionally, this section aims to demonstrate the

utility of GMC (which due to its ultra-efficient computation time, is more advantageous for use in a multiscale

framework than HFGMC or FEM) for progressive failure analysis. Although, the previous results are not

repeated on the following plots, those results are also considered when evaluating the performance of the

simpler models (whether HFGMC, GMC, or FEM).

A single fiber RUC, shown in Figure 31, representing a perfectly square packed architecture was modeled

using HFGMC and GMC. The 26 subcell x 26 subcell model contains 312 fiber subcells and 364 matrix

subcells. This configuration was chosen because it is one of the built-in architecture with in MAC/GMC,27

and since it is relatively coarse, is ideal for efficient multiscale simulations. Additionally, it should be noted

that due to the periodicity of the RUC, a single fiber RUC is identical to an RUC containing multiple square

packed fibers.

The results from HFGMC and GMC are compared to a more refined FEM model of the same fiber

architecture with the mesh shown in Figure 32. The FEM mesh consisted of 17,132 2D, quadrilateral,

generalized plane strain, reduced integration, CPEG4R elements, and periodic boundary conditions were

NASA/TM—2012-217649 13
applied to the edges of the RUC. This mesh was chosen to provide a fully-numerical baseline to compare the

results of the semi-analytical methods to.

The same five previous load cases were applied to all three models (HFGMC, GMC, and FEM). Additionally, the same elastic and fracture properties, as well as fiber volume fraction, as used in Section III were

also used in the single fiber simulations.

IV.B. Results for Square Packing Architecture

The transverse stress versus transverse strain response of the HFGMC, GMC, and FEM models of a single

fiber loaded with global transverse tensile strain is presented in Figure 33. The stiffnesses of all three models

correspond well, indicating that the previously observed discrepancy in the GMC stiffness is a facet of the

random architecture. Thus, due to the lack of shear coupling within the GMC formulation, GMC should

not be used to simulate disordered RUCs. Ultimate transverse stresses of 56.8 MPa, 54.6 MPa, and 51.3

MPa were predicted by the HFGMC, GMC, and FEM models, respectively. The ultimate stresses are higher

than what was obtained from the randomly arranged architectures. With the random arrangement of fibers,

two fibers will inevitably lie closer to each other than the fiber spacing in a square packed architecture.

Thus, the stress concentration between the more closely pack fibers will be higher, yielding a site for failure

initiation, followed by localization and subsequent propagation, at a lower global stress. In addition, the

FEM model was able to realize the stress concentration in the matrix better than HFGMC or GMC because

of the large number of matrix elements between the fiber boundary and RUC boundary. With HFGMC,

local fields are dependent on the subcell mesh (although not monotonically), much in the same way as

FEM solutions (this mesh dependency is not to be confused with pathological mesh dependence associated

with strain localization due to post-peak strain softening). Therefore, increasing of the number of subcells

between the fiber boundary and RUC boundary should lead) to a better resolution in the stress concentration

calculated using HFGMC. Alternatively and advantageously, the GMC solution is completely objective to

any refinements in the subcell mesh as long as the geometry of the RUC remains unchanged. Figure 34

displays the final mode I failure patterns for HFGMC, GMC and FEM solutions. All three models predicted

very similar failure modes.

In Figure 35 the transverse stress-strain response for a single fiber RUC subjected to global transverse

compressive strain is displayed. The non-linear behavior prior to sudden drops in the global stress obtained

from HFGMC closely match that from FEM. HFGMC predicted a severe stress drop at σ22 = -154.8 MPa

and 222 = -0.009, while FEM predicted a less severe stress drop at σ22 = -152.7 and 222 = -0.0092. The more

severe drop in the HFGMC stress state can be attributed to the presence of only one subcell between the

fiber boundary and the boundary of the RUC. GMC showed a similar trend as HFGMC and FEM, but the

response was much softer and the ultimate stress was -131 MPa. The transverse compressive failure modes

predicted by all three models, shown in Figure 36, are very similar. All three models predict a large amount

of mode II failure (due to local principle compression) at the fiber matrix interface on either side of the fiber.

These mode II regions are bridged by mode I crack bands on the top and/or bottom of the fiber. This failure

mode is consistent with some of the local failure surrounding individual fibers in Figure 13 and 14.

The transverse shear response of the single fiber RUCs are exhibited in Figure 37. There is far better

agreement between the stiffness of the GMC RUC and the HFGMC/FEM RUCs. However, the ultimate

shear stress predicted by GMC is higher than HFGMC and FEM: 49.2 MPa as compared to 40.8 MPa and

41.8 MPa. Moreover, the GMC RUC exhibits linear behavior up to failure, while HFGMC and FEM display

some nonlinearity. There is good agreement between the pre-peak response of the HFGMC and FEM RUC,

but the post-peak response of the FEM model is tougher. The failure patterns predicted with the three

models, which are displayed in Figure 38, show some variance. HFGMC predicts some distributed mode I

failure in the subcells at the top and bottom boundaries of the RUC accompanied by more substantial mode

II failure as well as a small degree of mode II failure in the matrix subcells between the closest of adjacent

fiber boundaries. In the GMC model complete matrix failure occurs in the matrix subcells adjacent to four

outermost fiber subcells, and the is a progression of mode II crack bands propagating vertically at the left

and right boundaries of the RUC. The FEM model predicts a moderately developed mode I crack bands in

the matrix at the fiber-matrix interface. A large amount of mode II failure can be seen in the FEM model

in an ‘L’-shaped pattern at the boundaries of the RUC. These failure modes are inconsistent with what was

predicted by the multiple fiber case, which exhibited no mode II crack bands and a mode I crack band that

propagated at an angle throughout the RUC.

NASA/TM—2012-217649 14
transverse tensile strain and transverse shear strain loading. There was some difference in the predicted peak

stresses between the models which is related to the variation of the predicted failure modes, given in Figure

40. HFGMC predicted mode I crack bands propagating around the fiber at the fiber-matrix interface, while

GMC predicted vertical mode II failure at the RUC boundaries, and the FEM model exhibited a single mode

I crack band at the fiber-matrix interface and into the matrix rich region near the bottom-left corner of the

RUC with some slight mode II failure throughout the RUC.

For the combined transverse compression-transverse shear load case, there is reasonable agreement between all three models (Figure 41). The ultimate shear stress calculated using FEM was slightly higher

than HFGMC and GMC, and the there was some disagreement in the post-peak toughnesses with the FEM

model exhibiting the toughest behavior and GMC the most brittle. The HFGMC and GMC failure patterns

where consistent, indicating mode II crack failure at the fiber-matrix interface with some mode I bridging,

while the FEM results showed some mode II failure propagating off the fiber-matrix interface at an angle

with some mode I crack bands at the fiber-matrix interfaces located at the top and bottom of the fiber.

V. Conclusions

In this paper, semi-analytical, mesh objective micromechanics methods, that are computationally efficient

(much faster that full finite element models) have been presented and validated against FEM based results

and when available, against experimental data. These methods (HFGMC and GMC) can be used for efficient

multi-scale, high-fidelity, progressive damage and failure (PDF) predictions of fiber reinforced composite

structures. Lack of mesh objectivity in damage and failure prediction methods, at all length scales, render

progressive damage and failure (PDF) tools to be “simulators” rather than“predictors”. Resolving stresses

and strains at the fiber-matrix level accurately and capturing their accumulated effect for upscaling to

the macroscale can be achieved at significantly reduced computational cost by using GMC/HFGMC in a

multiscale FEM based framework for PDF predictions.

The use of the smeared crack band theory within the HFGMC micromechanics model was verified for

transverse tension and compression against experimental data using an RUC containing 13 randomly arranged fiber. Furthermore, the implementation was validated against an analogous FEM model for transverse tension and compression, as well as transverse shear and combine transverse normal-transverse shear

loading cases. Good agreement was achieved between the HFGMC, FEM, and experimental results for these

load cases. A highly disordered RUC surpasses the intended applicability of GMC. Thus, as expected, GMC

performed poorly in predicting the stiffness and nonlinear behavior of the RUC, compared with the other

methods, as the shortcomings emerging from the lack of shear coupling (resulting from first order displacement approximation and displacement continuity restrictions) are further amplified when the packing of the

RUC is disordered. Additionally, the continuity conditions used to formulate the strain concentration matrix

in GMC creates a nonlocal relationship between the subcells in the RUC, and hence, GMC is not expected

to capture strain localization in a complex RUC well. The simulations were executed to determine a limit on

the error of the GMC solutions given an extreme case. The predicted failure modes obtained form HFGMC

and FEM corresponded well and represented physical failure modes observed from experiment.

A simpler, single fiber RUC (representing square packing architecture) was used to determine if accordance between GMC and the HFGMC/FEM could be achieved. In these models, the stiffness predictions

for GMC were greatly improved (compared with HFGMC/FEM). Additionally, there was an improvement

in the nonlinear results obtained from GMC. Global failure initiation was predicted reasonably. However,

GMC exhibited some error in post-peak failure progression. There was also some slight discrepancy between

the HFGMC and FEM results. In these analyses, a very refined FEM mesh was used to provide a trusted

baseline for the semi-analytical methods. Conversely, the HFGMC/GMC subcell meshes were relatively

coarse. It is expected that if the number of subcells between the fiber boundary and the periodic boundary

in the HFGMC model is increased, the error between the HFGMC and FEM results will be minimal.

Future work should study the effects of the size of the RUC containing randomly oriented fibers. Additionally, statistical variations in the fracture properties of the matrix should be used to introduce the

“randomness” of the multiple fiber RUC to a single fiber, square packed or dual-fiber hexagonally-packed

RUC. Moreover, the two cases presented here represent two extremes in disorder (or order). To clearly

distinguish the limitations of GMC as micromechanics model for predicting progressive failure, several intermediate levels of geometric perturbations within the RUC should be investigated. Finally, a methodology

to include the 3D effects, from planar loading of a composite lamina, in a more efficient 2D RUC should be

NASA/TM—2012-217649 15

Figure 39 contains plots of the global σ22-222 and τ23-γ23 responses of the RUCs subjected to combined
examined. By accomplishing these goals, a simple, efficient micromechanics methodology can be actualized

that contains all the physics of more encompassing models and will be readily suitable for computationally

efficient multi-scale applications.

NASA/TM—2012-217649 16
Property Value

Ef (GPa) 74.000

νf 0.200

Em (GPa) 4.65

νm 0.350

Gm 1.4

2

(βγ)

C

0.0143

γ

(βγ)

C

0.0237

G

(βγ)

IC (N/mm) 0.000563

G

(βγ)

IC (N/mm) 0.00385

Table 1: Elastic and fracture properties used in RUC simulation of 2-3 plane of a Silenka EGlass/MY750/HY917/DY063 composite lamina.

NASA/TM—2012-217649 17
Figure 1: Effects of mesh refinement on overall load-displacement response for a material exhibiting postpeak softening. Dashed lines indicate non-physical snapback.

Figure 2: Crack band domain Ω0 of width wc oriented normal to vector n within a continuum Ω.

NASA/TM—2012-217649 18
(a) (b)

Figure 3: Stress-strain (a) and traction-separation (b) laws governing material behavior. Failure energy

density WF is related to fracture toughness GC through the characteristic length.

Figure 4: Crack band embedded in discretized continuum. Magnified subcell displays crack band orientation

within subcell as well as characteristic length of subcell.

(a) Multiple transverse cracks.45 (b) Magnification of a single transverse crack.46

Figure 5: Transverse cracks in glass/epoxy laminates.

NASA/TM—2012-217649 19
(a) Macroscopic failure.47 (b) Localized shear bands and cracks.17

Figure 6: Transverse compressive failure of carbon fiber/epoxy laminates. Macroscopic failure path observed

at an angle. Magnification illustrates shear banding, cracking and fiber-matrix debonding.

Figure 7: HFGMC mesh used in simulations of RUC in 2-3 plane of lamina containing 13 fibers and a

fiber volume fraction of 58%. Subcells/elements occupied by fibers are colored blue, and subcells/elements

occupied by matrix are colored green.

NASA/TM—2012-217649 20
Figure 8: Transverse tensile stress versus transverse tensile strain of E-glass/Epoxy composite from experiment compared to HFGMC, GMC, and FEM micromechanics models containing 13 randomly placed

fibers.

(a) At ultimate stress (222 = 0.00248,

σ22 = 40.2 MPa).

(b) At second stress peak (222 =

0.00272, σ22 = 34.5 MPa).

(c) After load carrying capability of

has diminished (222 = 0.00288, σ22

= 2.33 MPa).

Figure 9: Mode I damage in RUC containing 13 randomly placed fibers subjected to applied transverse

tensile strain obtained using HFGMC.

NASA/TM—2012-217649 21
(a) At ultimate stress (222 =

0.00253, σ22 = 41.7 MPa).

(b) At second stress peak (222

= 0.0028, σ22 = 35.8 MPa).

(c) At third stress peak (222

= 0.00324, σ22 = 36.3 MPa).

(d) After load carrying capability of has diminished (222

= 0.0033, σ22 = 1.91 MPa).

Figure 10: Mode I damage in RUC containing 13 randomly placed fibers subjected to applied transverse

tensile strain obtained using FEM.

NASA/TM—2012-217649 22
(a) At first stress peak (222 = 0.0032,

σ22 = 41.5 MPa).

(b) At ultimate stress (222 = 0.0044,

σ22 = 46.7 MPa).

(c) At seventh stress peak (222 =

0.0057, σ22 = 32.1 MPa).

(d) After load carrying capability of

has diminished (222 = 0.01, σ22 =

3.64 MPa).

Figure 11: Mode I damage in RUC containing 13 randomly placed fibers subjected to applied transverse

tensile strain obtained using GMC.

NASA/TM—2012-217649 23
Figure 12: Transverse compressive stress versus transverse compressive strain of E-glass/Epoxy composite

from experiment compared to HFGMC, GMC, and FEM micromechanics models containing 13 randomly

placed fibers including mode I and mode II damage.

NASA/TM—2012-217649 24
(a) Mode II damage prior to ultimate

stress (222 = -0.0056, σ22 = -89.2

MPa).

(b) Mode II damage at ultimate stress

(222 = -0.00713, σ22 = -99.8 MPa).

(c) Mode II damage after load carrying capability has diminished (222 =

-0.0125, σ22 = -14.8 MPa).

(d) Mode I damage prior to ultimate

stress (222 = -0.0056, σ22 = -89.2

MPa).

(e) Mode I damage at ultimate stress

(222 = -0.00713, σ22 = -99.8 MPa).

(f) Mode I damage after load carrying capability has diminished (222 =

-0.0125, σ22 = -14.8 MPa).

(g) Mode I and mode II damage prior

to ultimate stress (222 = -0.0056, σ22

= -89.2 MPa).

(h) Mode I and mode II damage at

ultimate stress (222 = -0.00713, σ22

= -99.8 MPa).

(i) Mode I and mode II damage after load carrying capability has diminished (222 = -0.0125, σ22 = -14.8

MPa).

Figure 13: Mode I and Mode II damage in RUC containing 13 randomly placed fibers subjected to applied

transverse compressive strain obtained using HFGMC.

NASA/TM—2012-217649 25
(a) Mode II damage prior

to ultimate stress (222 = -

0.00582, σ22 = -94.1 MPa).

(b) Mode II damage at ultimate stress (222 = -0.00755,

σ22 = -108.3 MPa).

(c) Mode I damage prior

to ultimate stress (222 = -

0.00582, σ22 = -94.1 MPa).

(d) Mode I damage at ultimate stress (222 = -0.00755,

σ22 = -108.3 MPa).

(e) Mode I and mode II damage prior to ultimate stress

(222 = -0.00582, σ22 = -94.1

MPa).

(f) Mode I and mode II damage at ultimate stress (222 =

-0.00755, σ22 = -108.3 MPa).

Figure 14: Mode I and Mode II damage in RUC containing 13 randomly placed fibers subjected to applied

transverse compressive strain obtained using FEM.

NASA/TM—2012-217649 26
(a) Mode II damage prior to ultimate

stress (222 = -0.00429, σ22 = -55.4

MPa).

(b) Mode II damage at ultimate stress

(222 = -0.00528, σ22 = -62.4 MPa).

(c) Mode II damage after load carrying capability has diminished (222 =

-0.0125, σ22 = -0.65 MPa).

Figure 15: Mode II damage in RUC containing 13 randomly placed fibers subjected to applied transverse

compressive strain obtained using GMC.

Figure 16: Transverse compressive stress versus transverse compressive strain of E-glass/Epoxy composite

from experiment compared to HFGMC, GMC, and FEM micromechanics models containing 13 randomly

placed fibers with mode I damage disabled.

NASA/TM—2012-217649 27
(a) 222 = -0.0056, σ22 = -90.4 MPa. (b) 222 = -0.012, σ22 = -145.8 MPa. (c) 222 = -0.014, σ22 = -153.5 MPa.

Figure 17: Mode II damage in RUC containing 13 randomly placed fibers subjected to applied transverse

compressive strain obtained using HFGMC with mode I damage disabled.

(a) 222 = -0.00564, σ22 = -

91.1 MPa.

(b) 222 = -0.012, σ22 = -156.0

MPa.

(c) 222 = -0.014, σ22 = -167.5

MPa.

Figure 18: Mode I damage in RUC containing 13 randomly placed fibers subjected to applied transverse

compressive strain obtained using FEM with mode I damage disabled.

NASA/TM—2012-217649 28
Figure 19: Transverse shear stress versus transverse shear strain of E-glass/Epoxy composite from HFGMC,

GMC, and FEM micromechanics models containing 13 randomly placed fibers.

NASA/TM—2012-217649 29
(a) At initiation (γ23 = 0.0043, τ23 =

26.2 MPa).

(b) At ultimate stress (γ23 = 0.0053,

τ23 = 33.4 MPa).

(c) After RUC has lost all load carrying capability (γ23 = 0.0064, τ23 =

4.6 MPa).

Figure 20: Mode I damage in RUC containing 13 randomly placed fibers subjected to applied transverse

shear strain obtained using HFGMC.

(a) At first stress peak (γ23

= 0.00431, τ23 = 28.1 MPa).

(b) At ultimate stress (γ23 =

0.00543, τ23 = 33.3 MPa).

(c) After load carry capability

of RUC has diminished (γ23

= 0.00544, τ23 = 16.8 MPa).

Figure 21: Mode I damage in RUC containing 13 randomly placed fibers subjected to applied transverse

shear strain obtained using FEM.

Figure 22: Mode I damage in RUC containing 13 randomly placed fibers subjected to applied transverse

shear strain obtained using GMC at, and subsequent to ultimate load.

NASA/TM—2012-217649 30
(a) Transverse stress versus transverse strain. (b) Transverse shear stress versus transverse shear strain.

Figure 23: Response of E-glass/Epoxy composite RUC subjected to applied transverse shear strain and

transverse tension calculated from HFGMC, GMC, and FEM micromechanics models containing 13 randomly

placed fibers.

(a) Just after ultimate stress (222

= 0.00098, σ22 = 7.55 MPa γ23 =

0.00488, τ23 = 23.6 MPa).

(b) Just after second shear stress

peak (222 = 0.0013, σ22 = -5.75 MPa

γ23 = 0.0065, τ23 = 20.5 MPa).

(c) Final (222 = 0.0025, σ22 = -13.8

MPa γ23 = 0.0125, τ23 = 11.7 MPa).

Figure 24: Mode I damage in RUC containing 13 randomly placed fibers subjected to applied transverse

tensile strain and transverse shear strain obtained using HFGMC.

NASA/TM—2012-217649 31
(a) At first stress peak (222 =

0.00084, σ22 = 13.4 MPa γ23

= 0.0042, τ23 = 27.0 MPa).

(b) At ultimate stress (222 =

0.00096, σ22 = 13. MPa γ23

= 0.0048, τ23 = 27.1 MPa).

(c) Final (222 = 0.0016, σ22 =

-3.89 MPa γ23 = 0.008, τ23 =

3.25 MPa).

Figure 25: Mode I damage in RUC containing 13 randomly placed fibers subjected to applied transverse

tensile strain and transverse shear strain obtained using FEM.

(a) Mode I damage. (b) Mode II damage.

Figure 26: Damage in RUC containing 13 randomly placed fibers subjected to applied transverse tensile

strain and transverse shear strain obtained using GMC at end of simulation.

(a) Transverse stress versus transverse strain. (b) Transverse shear stress versus transverse shear strain.

Figure 27: Response of E-glass/Epoxy composite RUC containing 13 randomly placed fibers subjected

to applied transverse shear strain and transverse compression calculated from HFGMC, GMC, and FEM

micromechanics models.

NASA/TM—2012-217649 32
(a) Mode I damage at first transverse

stress peak (222 = -0.00405, σ22 = -

64.7 MPa, γ23 = 0.0054, τ23 = 31.1

MPa).

(b) Mode I damage at ultimate stress

(222 = -0.0069, σ22 = -77.3 MPa, γ23

= 0.0092, τ23 = 25.5 MPa).

(c) Mode I damage after load carrying capability has diminished (222 =

-0.015, σ22 = -10.1 MPa, γ23 = 0.02,

τ23 = 5.3 MPa).

(d) Mode II damage at first transverse

stress peak (222 = -0.00405, σ22 = -

64.7 MPa, γ23 = 0.0054, τ23 = 31.1

MPa).

(e) Mode II damage at ultimate stress

(222 = -0.0069, σ22 = -77.3 MPa, γ23

= 0.0092, τ23 = 25.5 MPa).

(f) Mode II damage after load carrying capability has diminished (222 =

-0.015, σ22 = -10.1 MPa, γ23 = 0.02,

τ23 = 5.3 MPa).

(g) Mode I and mode II at first transverse stress peak (222 = -0.00405, σ22

= -64.7 MPa, γ23 = 0.0054, τ23 = 31.1

MPa).

(h) Mode I and mode II damage at

ultimate stress (222 = -0.0069, σ22 =

-77.3 MPa, γ23 = 0.0092, τ23 = 25.5

MPa).

(i) Mode I and mode II damage after load carrying capability has diminished (222 = -0.015, σ22 = -10.1

MPa, γ23 = 0.02, τ23 = 5.3 MPa).

Figure 28: Mode I and Mode II damage in RUC containing 13 randomly placed fibers subjected to applied

transverse compressive strain and transverse shear strain obtained using HFGMC.

NASA/TM—2012-217649 33
(a) Mode I damage at first

stress peak (222 = -0.0038,

σ22 = -61.6 MPa γ23 =

0.00507, τ23 = 29.8 MPa).

(b) Mode I damage at end

of simulation (222 = -0.0052,

σ22 = -73.3 MPa γ23 =

0.00696, τ23 = 29.2 MPa).

(c) Mode II damage at first

stress peak (222 = -0.0038,

σ22 = -61.6 MPa γ23 =

0.00507, τ23 = 29.8 MPa).

(d) Mode II damage at

end of simulation (222 = -

0.0052, σ22 = -73.3 MPa γ23

= 0.00696, τ23 = 29.2 MPa).

(e) Mode I and mode II damage at first stress peak (222 =

-0.0038, σ22 = -61.6 MPa γ23

= 0.00507, τ23 = 29.8 MPa).

(f) Mode I and mode II damage at end of simulation (222

= -0.0052, σ22 = -73.3 MPa

γ23 = 0.00696, τ23 = 29.2

MPa).

Figure 29: Mode I and mode II damage in RUC containing 13 randomly placed fibers subjected to applied

transverse tensile strain and transverse shear strain obtained using FEM.

NASA/TM—2012-217649 34
(a) At ultimate stress (222 = -0.00435,

σ22 = -53.7 MPa γ23 = 0.0058, τ23 =

22.8 MPa).

(b) At second shear stress peak (222

= -0.0048, σ22 = -48.1 MPa γ23 =

0.0064, τ23 = 22.7 MPa).

(c) Final (222 = -0.015, σ22 = -1.95

MPa γ23 = 0.02, τ23 = 0.83 MPa).

Figure 30: Mode II damage in RUC containing 13 randomly placed fibers subjected to applied transverse

tensile strain and transverse shear strain obtained using GMC.

Figure 31: HFGMC and GMC RUC subcell mesh used for single fiber, square packing architecture simulations.

Figure 32: FEM RUC subcell mesh used for single fiber, square packing architecture simulations.

NASA/TM—2012-217649 35
Figure 33: Transverse tensile stress versus transverse tensile strain of E-glass/Epoxy composite RUC containing a single fiber from HFGMC, GMC, and FEM micromechanics models.

(a) HFGMC. (b) GMC. (c) FEM.

Figure 34: Final mode I damage in single fiber RUC subjected to applied transverse tensile strain.

NASA/TM—2012-217649 36
Figure 35: Transverse compressive stress versus transverse compressive strain of E-glass/Epoxy composite

RUC containing a single fiber from HFGMC, GMC, and FEM micromechanics models.

NASA/TM—2012-217649 37
(a) Mode I, HFGMC. (b) Mode I, GMC. (c) Mode I, FEM.

(d) Mode II, HFGMC. (e) Mode II, GMC. (f) Mode II, FEM.

(g) Mode I and mode II, HFGMC. (h) Mode I and mode II, GMC. (i) Mode I and mode II,

FEM.

Figure 36: Final damage in single fiber RUC subjected to applied transverse compressive strain.

NASA/TM—2012-217649 38
Figure 37: Transverse shear stress versus transverse shear strain of E-glass/Epoxy composite RUC containing

a single fiber from HFGMC, GMC, and FEM micromechanics models.

NASA/TM—2012-217649 39
(a) Mode I, HFGMC. (b) Mode I, GMC. (c) Mode I, FEM.

(d) Mode II, HFGMC. (e) Mode II, GMC. (f) Mode II, FEM.

(g) Mode I and mode II, HFGMC. (h) Mode I and mode II, GMC. (i) Mode I and mode II,

FEM.

Figure 38: Final damage in single fiber RUC subjected to applied transverse shear strain.

NASA/TM—2012-217649 40
(a) Transverse stress versus transverse strain. (b) Transverse shear stress versus transverse shear strain.

Figure 39: Response of E-glass/Epoxy composite RUC containing a single fiber subjected to applied transverse shear strain and transverse tension calculated from HFGMC, GMC, and FEM micromechanics models.

NASA/TM—2012-217649 41
(a) Mode I, HFGMC. (b) Mode I, GMC. (c) Mode I, FEM.

(d) Mode II, HFGMC. (e) Mode II, GMC. (f) Mode II, FEM.

(g) Mode I and mode II, HFGMC. (h) Mode I and mode II, GMC. (i) Mode I and mode II,

FEM.

Figure 40: Final damage in single fiber RUC subjected to applied transverse shear strain and transverse

tension.

NASA/TM—2012-217649 42
(a) Transverse stress versus transverse strain. (b) Transverse shear stress versus transverse shear strain.

Figure 41: Response of E-glass/Epoxy composite RUC containing a single fiber subjected to applied transverse shear strain and transverse compression calculated from HFGMC, GMC, and FEM micromechanics

models.

NASA/TM—2012-217649 43
(a) Mode I, HFGMC. (b) Mode I, GMC. (c) Mode I, FEM.

(d) Mode II, HFGMC. (e) Mode II, GMC. (f) Mode II, FEM.

(g) Mode I and mode II, HFGMC. (h) Mode I and mode II, GMC. (i) Mode I and mode II,

FEM.

Figure 42: Final damage in single fiber RUC subjected to applied transverse shear strain and transverse

compression.

NASA/TM—2012-217649 44
References

1Aboudi, J., Arnold, S. M., and Bednarcyk, B. A., Micromechanics of Composite Materials: A Generalized Multiscale

Analysis Approach, Elsevier, 2012, in press.

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

3Paley, M. and Aboudi, J., “Micromechanical analysis of composites by the generalized cells model,” Mechanics of Materials, Vol. 14, 1992, pp. 127–139.

4Aboudi, J., “Micormechanical analysis of thermo-inelastic multiphase short-fiber composites,” Compos. Eng., Vol. 5,

No. 7, 1995, pp. 839–850.

5Aboudi, J., Pindera, M.-J., and Arnold, S. M., “Linear thermoelastic higher-order theory for periodic multiphase materials,” J. Appl. Mech., Vol. 68, 2001, pp. 697–707.

6Haj-Ali, R. and Aboudi, J., “Nonlinear micromechanical formulation of the high fidelity generalized method of cells,”

Int. J. Solids Struct., Vol. 46, 2009, pp. 2577–2592.

7Aboudi, J., Pindera, M.-J., and Arnold, S. M., “Higher-order theory for periodic multiphase materials with inelastic

phases,” Int. J. Plast., Vol. 19, 2003, pp. 805–847.

8Bednarcyk, B. A., Arnold, S. M., Aboudi, J., and Pindera, M.-J., “Local field effects in titanium matrix composites

subject to fiber-matrix debonding,” Int. J. Plast., Vol. 20, 2004, pp. 1707–1737.

9Bednarcyk, B. A., Aboudi, J., and Arnold, S. M., “Micromechanics modeling of composites subjected to multiaxial

progressive damage in the constituents,” AIAA J., Vol. 48, 2010, pp. 1367–1378.

10Pindera, M. J. and Bednarcyk, B. A., “An efficient implementation of the GMC micromechanics model for multi-phased

materials with complex microstructures,” NASA/CR 1997-202350, 1997.

11Baˇzant, Z. P. and Oh, B. H., “Crack band theory for fracture of concrete,” Mater. and Struct., Vol. 16, 1983, pp. 155–77.

12de Borst, R. and Nauta, P., “Non-orthogonal cracks in a smeared finite element model,” Eng. Comput., Vol. 2, 1985,

pp. 35–46.

13Rots, J. G. and de Borst, R., “Analysis of mixed-mode fracture in concrete,” J. Eng. Mech., Vol. 113, No. 11, 1987,

pp. 1739–1758.

14Camanho, P. P., Maim ́ı, P., and D ́avila, C. G., “Prediction of size effects in notched laminates using continuum damage

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

15Spencer, B. W., Finite Elements with Embedded Discontinuities for Modeling Reinforced Concrete Members, Ph.D.

thesis, Brigham Young University, Provo, UT, 2002.

16Pineda, E. J. and Waas, A. M., “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, 2011.

17Gonzal ́ez, C. and Llorca, J., “Mechanical behavior of unidirectional fiber-reinforced polymers under transverse compression: Microscopic mechanisms and modeling,” Compos. Sci. Technol., Vol. 67, 2007, pp. 2795–2806.

18Totry, E., Gonzal ́ez, C., and Llorca, J., “Failure locus of fiber-reinforced composites under transverse compression and

out-of-plane shear,” Compos. Sci. and Technol., Vol. 68, 2008, pp. 829–839.

19Totry, E., Gonzalez, C., and Llorca, J., “Effect of fiber, matrix and interface properties on the in-plane shear deformation

of carbon-fiber reinforced composites,” Compos. Sci. and Technol., Vol. 70, 2010, pp. 970–980.

20Pineda, E. J., Bednarcyk, B. A., Waas, A. M., and Arnold, S. M., “Implementation of a smeared crack band model in a

micromechanics framework,” NASA/TM 2012-217603, 2012.

21Baˇzant, Z. and Cedolin, L., “Blunt crack band propagation in finite element analysis,” J. Eng. Mech. Div.-ASCE,

Vol. 105, 1979, pp. 297–315.

22Pietruszczak, S. and Mroz, Z., “Finite element analysis of deformation of strain-softening materials,” Int. J. Numer.

Methods Eng., Vol. 17, 1981, pp. 327–334.

23deBorst, R., “Computation of post-bifurcation and post-failure behavior of strain-softening solids,” Comput. Struct.,

Vol. 25, 1987, pp. 211–224.

24Baˇzant, Z. P. and Cedolin, L., Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories, Oxford University Press, New York, Oxford, 1991.

25Baˇzant, Z. P., “Crack band model for fracture of geomaterials,” Proceedings of the 4th International Conference on

Numerical Methods in Geomechanics, Edmonton, Canada, 1982.

26Baˇzant, Z. P. and Cedolin, L., “Finite element modeling of crack band propagation,” J. of Struc. Eng., Vol. 109, 1983,

pp. 69–92.

27Bednarcyk, B. A. and Arnold, S. M., “MAC/GMC 4.0 User’s Manual - Keywords Manual,” Nasa tm-2002-212077/vol2,

2002.

28Bednarcyk, B. A. and Arnold, S. M., “MAC/GMC 4.0 User’s Manual - Example Problems Manual,” Nasa tm-2002-

212077/vol3, 2002.

29Tay, T. E., Liu, G., C., V. B., Tan, Sun, X. S., and Pham, D. C., “Progressive failure analysis of composites,” J. Compos.

Mater., Vol. 42, 2008, pp. 1921–1966.

30Hoek, E. and Bieniawski, Z. T., “Brittle rock fracture propagation in rock under compression,” Int. J. Fract., Vol. 1,

No. 3, 1965, pp. 137–155.

31Horii, H. and Nemat-Nasser, S., “Brittle failure in compression: Splitting, faulting and brittle-ductile transition,” Philos.

T. R. Soc. S.-A., Vol. 319, No. 1549, 1986, pp. 337–374.

32Ashby, M. F. and Sammis, C. G., “The damage mechanics of brittle solids in compression,” Pure Appl. Geophys., Vol. 133,

No. 3, 1990, pp. 489–520.

33Chen, W. and Ravichandran, G., “Failure mode transition in creamics under dynamic multiaxial compression,” Int. J.

Fract., Vol. 101, 2000, pp. 141–159.

NASA/TM—2012-217649 45
34Taliercio, A. and Sagramoso, P., “Uniaxial strength of polymeric-matrix fibrous composites predicted through a homogenization approach,” Int. J. Solids Struct., Vol. 32, No. 14, 1995, pp. 2095–2123.

35Puck, A. and Sch ̈urmann, H., “Failure analysis of FRP laminates by means of physically based phenomenological models,”

Comps. Sci. Technol., Vol. 58, 1998, pp. 1045–1067.

36Puck, A. and Sch ̈urmann, H., “Failure analysis of FRP laminates by means of physically based phenomenological models,”

Comps. Sci. Technol., Vol. 62, 2002, pp. 1633–1622.

37Pinho, S. T., D ́avila, C. G., Camanho, P. P., and Iannucci, L., “Failure models and criteria for FRP under in-plane or

three-dimensional stress states including shear non-linearity,” NASA/TM 2005-213530, 2005.

38Haj-Ali, R. and Aboudi, J., “Formulation of the high-fidelity generalized method of cells with arbitrary cell geometry for

refined micromechanics and damage in composites,” Int. J. Solids Struct., Vol. 47, 2010, pp. 3447–3461.

39Pineda, E. J., Waas, A. M., Bednarcyk, B. A., and Collier, C. S., “An efficient semi-analytical framework for micromechanical modeling of transverse cracks in fiber-reinforced composites,” 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural

Dynamics, and Materials Conference, 12-15 April 2010.

40Pineda, E. J., Waas, A. M., Bednarcyk, B. A., and Collier, C. S., “Multiscale Modeling of Progressive Microdamage

and Transverse Cracking in Fiber Reinforced Laminates,” Presented at the 16th US National Conference of Theoretical and

Applied Mechanics, 27 June - 2 July 2010.

41Hinton, M. J., Kaddour, A. S., and Soden, P. D., editors, Elsevier, New York, 2004.

42Arnold, S. M., Pindera, M. J., and Wilt, T. E., “Influence of fiber architecture on the inelastic response of metal matrix

composites,” Int. J. Plasticity, Vol. 12, No. 4, 1996, pp. 507–545.

43Ng, 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., Vol. 70, No. 7, 2010, pp. 1126–1134.

44Abaqus, Abaqus User’s Manual, Vol. 1-3, Version 6.10-1 , Dassault Syst`emes Simulia Corp., Providence, RI, 2008.

45Roberts, S. J., Modelling of Microcracking in Composite Materials, Ph.D. thesis, The University of Newcastle upon Tyne,

Newcastle upon Tyne, United Kingdom, 2000.

46Gamstedt, E. K. and Sj ̈ogren, B. A., “Micromechanisms in tension-compression fatigue of composite laminates containing

transverse plies,” Compos. Sci. Technol., Vol. 59, 1999, pp. 167–178.

47Aragon ́es, D., Fracture Micromechanisms in C/Epoxy Composites under Transverse Compression, Ph.D. thesis, Universidad Polit ́ecnica de Madrid, Madrid, Spain, 2007.

NASA/TM—2012-217649 46

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

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Progressive Failure of a Unidirectional Fiber-Reinforced Composite Using the Method of 

Cells: Discretization Objective Computational Results 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Pineda, Evan, J.; Bednarcyk, Brett, A.; Waas, Anthony, M.; Arnold, Steven, M. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 984754.02.07.03.16.03.02 

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

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

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

The smeared crack band theory is implemented within the generalized method of cells and high-fidelity generalized method of cells 

micromechanics models to capture progressive failure within the constituents of a composite material while retaining objectivity with 

respect to the size of the discretization elements used in the model. An repeating unit cell containing 13 randomly arranged fibers is modeled 

and subjected to a combination of transverse tension/compression and transverse shear loading. The implementation is verified against 

experimental data (where available), and an equivalent finite element model utilizing the same implementation of the crack band theory. To 

evaluate the performance of the crack band theory within a repeating unit cell that is more amenable to a multiscale implementation, a single 

fiber is modeled with generalized method of cells and high-fidelity generalized method of cells using a relatively coarse subcell mesh which 

is subjected to the same loading scenarios as the multiple fiber repeating unit cell. The generalized method of cells and high-fidelity 

generalized method of cells models are validated against a very refined finite element model. 

15. SUBJECT TERMS

Micromechanics; 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

54 

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