Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Jacob Aboudi

Tel Aviv University, Ramat Aviv, Israel

Steven M. Arnold

Glenn Research Center, Cleveland, Ohio

Analysis of Fiber Clustering in Composite Materials 

Using High-Fidelity Multiscale Micromechanics

NASA/TM—2015-218838

July 2015
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated 

to the advancement of aeronautics and space science. 

The NASA Scientific and Technical Information (STI) 

Program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Officer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI Program provides access 

to the NASA Technical Report Server—Registered 

(NTRS Reg) and NASA Technical Report Server—

Public (NTRS) thus providing one of the largest 

collections of aeronautical and space science STI in 

the world. Results are published in both non-NASA 

channels and by NASA in the NASA STI Report 

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of significant 

scientific and technical data and information 

deemed to be of continuing reference value. 

NASA counter-part of peer-reviewed formal 

professional papers, but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., “quick-release” reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis.

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA.

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission.

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 757-864-6500

• Telephone the NASA STI Information Desk at

757-864-9658

• Write to:

NASA STI Program

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199
Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Jacob Aboudi

Tel Aviv University, Ramat Aviv, Israel

Steven M. Arnold

Glenn Research Center, Cleveland, Ohio

Analysis of Fiber Clustering in Composite Materials 

Using High-Fidelity Multiscale Micromechanics

NASA/TM—2015-218838

July 2015

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Acknowledgments

The authors gratefully acknowledge Drs. Bradley A. Lerch and Pappu L.N. Murthy of NASA Glenn Research Center

who provided the composite micrographs and the image pixel analysis results used herein.

Available from

Trade names and trademarks are used in this report for identification 

only. Their usage does not constitute an official endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

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

NASA STI Program

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22161

703-605-6000

This report is available in electronic form at http://www.sti.nasa.gov/ and http://ntrs.nasa.gov/
NASA/TM—2015-218838 1 

Analysis of Fiber Clustering in Composite Materials Using 

High-Fidelity Multiscale Micromechanics

Brett A. Bednarcyk

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135 

Jacob Aboudi 

Tel Aviv University

Ramat Aviv 69978, Israel 

Steven M. Arnold 

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135 

Abstract

A new multiscale micromechanical approach is developed for the prediction of the behavior of fiber 

reinforced composites in presence of fiber clustering. The developed method is based on a coupled twoscale implementation of the High-Fidelity Generalized Method of Cells theory, wherein both the local and 

global scales are represented using this micromechanical method. Concentration tensors and effective 

constitutive equations are established on both scales and linked to establish the required coupling, thus 

providing the local fields throughout the composite as well as the global properties and effective 

nonlinear response. Two nondimensional parameters, in conjunction with actual composite micrographs, 

are used to characterize the clustering of fibers in the composite. Based on the predicted local fields, 

initial yield and damage envelopes are generated for various clustering parameters for a polymer matrix 

composite with both carbon and glass fibers. Nonlinear epoxy matrix behavior is also considered, with 

results in the form of effective nonlinear response curves, with varying fiber clustering and for two sets of 

nonlinear matrix parameters.

Introduction

Clustering in composite materials refers to the aggregation of the constituents such that the composite 

does not exhibit uniform microstructural distribution. Composite microstructures are sensitive to the 

specific manufacturing processes utilized, and because the processing cannot be perfectly controlled, 

some degree of clustering will always be present. Clustering is more prevalent in composites with low 

fiber or inclusion volume fractions as there is more volume available in which one constituent may be 

segregated. In fact, in nanocomposites, whose volume fractions are typically very low, clustering is a 

dominant phenomenon, which strongly influences the properties and performance of the composite 

material (c.f., Shaffer and Windle (1999), Vigolo et al. (2000), Shi et al. (2004)). 

Several investigators have presented methods to characterize the degree of clustering in composite 

materials (c.f., Guild and Summerscales (1993), Pyrz (1994), Scanlon et. al. (2003), Al-Ostaz et al. 

(2007), Vaughan and McCarthy (2010), Wilding and Fulwood (2011), Zangenberg et al. (2013)). In terms 

of modeling the effects of fiber clustering in composites, most investigations have been finite element 

based, as this approach can explicitly account for the composite microstructure. Examples include 
NASA/TM—2015-218838 2 

Sorensen and Talreja (1993), Ghosh et al. (1997), Yang et al. (1997, 2000), Zeman and Sejnoha (2001), 

Wongsto and Li (2005), Melro et al. (2008), Abhilash et al. (2011), and Romanov et al. (2013). 

Additionally, the Generalized Method of Cells semi-analytical micromechanics theory has been employed 

to model random fiber distributions in continuous composites with a moving-window technique (Baxter 

and Graham (2000), Graham-Brady et al. (2003), Baxter et al. (2005), Acton and Graham-Brady (2010)). 

In the present investigation, the triply-periodic High-Fidelity Generalized Method of Cells (HFGMC) 

(Aboudi et al., 2013) is first enhanced to enable coupled multiscale analysis, wherein both the local and 

global scales are synergistically linked. The HFGMC method determines the strain (or stress) 

concentration tensors, which are used to establish the macroscopic constitutive equations of the 

composite, and also to provide the local stress and strain fields throughout the composite. This enables the 

prediction of not only effective composite properties, but also the effective nonlinear response (due to 

local inelasticity and damage) in response to full multiaxial loading, as well as generation of initial yield 

and damage envelopes. The present multiscale HFGMC implementation considers an arbitrary number of 

local repeating unit cells (RUCs), which consist of the composite constituent materials (e.g., fiber and 

matrix). These local RUCs are then arranged in a global RUC, whose behavior represents that of the 

multiscale composite material. The scales are explicitly coupled so that any local nonlinearity (which can 

be predicted using the local fields) is homogenized and passed to the global scale analysis. The examples 

presented here utilize a two scale analysis, but the methodology could be generalized to admit an arbitrary 

number of scales, as has been done in the case of the Multiscale Generalized Method of Cells (Liu et al., 

2011, Liu and Arnold, 2013; Bednarcyk et al., 2015). 

Herein, the multiscale HFGMC approach is employed to examine the effects of fiber clustering in 

continuous fiber reinforced composites. A key aspect of the approach taken involves segregating the 

composite microstructure into two zones, A and B, see Figure 1. Each zone consists of a composite 

material with its own fiber volume fraction, thus, by varying the fiber volume fractions and the 

arrangement/size of zones A and B, fiber clustering can be represented. For example, if the fiber volume 

fraction in zone A is higher than that in zone B, the fiber-rich zone A represents a clustered region.

Conversely, if the fiber volume fraction in zone B is higher than that in zone A, zone A will be a matrixrich region as compared to zone B. It should be noted that the present method admits an arbitrary number 

of zones. The alternative single scale HFGMC or finite element approach to modeling fiber clustering 

would involve full discretization of the actual composite microstructure. Obviously, this implies the use 

of a very dense mesh with great potential for meshing issues in the local fiber-matrix regions. In contrast, 

the proposed multiscale HFGMC methodology captures the effects of fiber clustering through local scale 

fiber volume fraction variations. Thus coarse meshes can be used on both scales, avoiding these meshing 

concerns while balancing efficiency and fidelity. Further, while mean field approaches (e.g., Mori-Tanaka 

method) can predict reliable effective properties, the predicted mean fields do not capture the variations of 

the stress fields in the composite constituents, even in the case of uniform fiber distribution. Utilization of 

mean fields to predict initial yield envelopes was shown to be problematic via comparison to finite 

element micromechanics analyses by Pindera and Aboudi (1988) in the case of uniform fiber distribution.

In the presence of fiber clustering, the matrix stress field variations are clearly important, thus motivating 

the use of a higher fidelity micromechanical approach such as HFGMC. It should be noted, however, that 

lower fidelity methods with better computational efficiency, such as the Generalized Method of Cells or 

the Mori-Tanaka Method, could be used to generate the results like those presented herein. The impact on 

these results of the better local fields given by HFGMC versus other methods has not been investigated. 
NASA/TM—2015-218838 3 

Figure 1.—Schematic drawing of a composite 

divided into two zones, A and B, which have 

different fiber volume fractions.

In the present investigation, actual micrographs from polymer matrix composites are used to 

determine these two zones within the composite RUC. The fiber clustering is characterized by two 

nondimensional parameters, one representing the volume fraction of zone A within the composite, and the 

other representing the fraction of the fibers that are located in zone A. Studies are performed by varying 

these parameters and observing their impact on the predicted initial global-scale yield and damage 

envelopes, as well as the nonlinear global stress-strain response. It should be noted that the availability of 

the concentration tensors in HFGMC enables straightforward and efficient generation of these initial 

envelopes. In contrast, the use of the standard finite element method would require repeated application of 

many loading combinations to trace out these surfaces.

The remainder of this paper is organized as follows. First, the multiscale HFGMC theory is presented 

for the generally triply-periodic case. Due to the possible nonlinearity of the constituent materials, an 

incremental (tangential) formulation is employed. The local and global concentration tensors and 

constitutive equations, as well as the coupling between them, are established. Next the methodology for 

determination of the initial yield and damage envelopes is described, followed by the method used to 

model fiber clustering. The Results and Discussion section details the processing of the composite 

micrographs to determine the analyzed composite RUCs. Results in the form of predicted effective 

properties, initial yield and damage envelopes, and macroscopic nonlinear stress-strain curves are 

presented for carbon/epoxy, glass/epoxy, and two different nonlinear representations of the epoxy matrix. 

The Conclusion section summarizes the paper and offers future research directions.

Multiscale High-Fidelity Generalized Method of Cells

In order to model the clustering of fibers in a composite material, the High-Fidelity Generalized 

Method of Cells (HFGMC) micromechanical model has been implemented in a two-scale framework.

These two scales are referred to as global and local. The global scale represents a repeating unit cell of a 

periodic composite material, whose constituents are themselves periodic composite materials. Thus, the 

local scale represents the RUCs present within the global scale constituents, see Figure 2. 
NASA/TM—2015-218838 4 

Figure 2.—(a) A multiphase composite with triply-periodic microstructure defined with respect to global coordinates 

(X1, X2, X3). (b) Its repeating unit cell (RUC) is represented with respect to the coordinates (Y1, Y2, Y3) and is 

divided into NΑ, NΒ, and NΓ global-scale subcells, in the Y1, Y2, and Y3 directions, respectively. (c) The materials 

within the global-scale subcells (ΑΒΓ) are represented by local-scale RUCs, defined with respect to local

coordinates (y1, y2, y3), and is divided into Nα, Nβ, and Nγ local-scale subcells. 

Global Scale Analysis

The HFGMC theory, which has been fully described by Aboudi et al. (2013), considers a composite 

material with triply-periodic microstructure, Figure 2(a), wherein periodicity conditions are enforced in 

all three Cartesian coordinate directions. The global repeating unit cell (RUC), Figure 2(b), defined with 

respect to local coordinates (Y1, Y2, Y3), is divided into NΑ, NΒ, and NΓ global subcells in the Y1, Y2, and Y3

directions, respectively. Each global subcell is labeled by the indices (ΑΒΓ) with Α = 1 ,...,NΑ, 

Β = 1,...,NΒ and Γ = 1,...,NΓ, and may contain a distinct homogeneous material or a composite material. 

The dimensions of the RUC are D, H, and L, whereas the dimensions of global subcell (ΑΒΓ) in the Y1, 

Y2, and Y3 directions are denoted by DΑ, HΒ, and LΓ, respectively. A coordinate system ),,(

)(

3

)(

2

)(

2

ΓΒΑ YYY

is introduced in each subcell whose origin is located at its center. The global subcell nonlinear elastic 

constitutive equation of the anisotropic material is given in an incremental form by, 

(a) 

(b) 

(c) 
NASA/TM—2015-218838 5 

ΑΒΓ ΑΒΓ ΑΒΓ)()()( ε∆=σ∆ ij ijkl kl C (1)

where ΑΒΓ)( σ∆ ij , ΑΒΓ)( ε∆ kl , and ΑΒΓ)(

ijkl C are the components of the stress increment, strain increment, and

instantaneous (tangent) stiffness tensors of global subcell (ΑΒΓ), respectively.

The basic assumption in HFGMC is that the increment of the displacement vector ΑΒΓ)( ∆Ui in each 

global subcell is represented as a second-order expansion in terms of its coordinates ),,(

)(

3

)(

2

)(

1

ΓΒΑ YYY , as 

follows, 

)(

)002(

2

2)(

3

)(

)020(

2

2)(

2

)(

)200(

2

2)(

1

)(

)001(

)(

3

)(

)010(

)(

2

)(

)100(

)(

1

)(

)000(

)(

4

3

2

1

4

3

2

1

4

3

2

1 Α Α ΑΒΓ Β Β ΑΒΓ Γ Γ ΑΒΓ

ΑΒΓ ΑΒΓ Α ΑΒΓ Β ΑΒΓ Γ ΑΒΓ

∆















+∆ − 















+∆ − 















+ −

∆+∆+∆+∆+ε∆=∆

i i i

i i i i i ij j

W

L

W Y

H

W Y

D

Y

U WYWYWYWX

(2)

where ε∆ ij are the applied (external) average strain increments, and the unknown terms )(

)(

ΑΒΓ ∆Wi lmn must 

be determined from the fulfillment of the equilibrium conditions, the periodic boundary conditions, and 

the interfacial continuity conditions of displacements and tractions between global subcells. The periodic 

boundary conditions ensure that the increments of displacement and traction at opposite surfaces of the 

global RUC are identical. A principal ingredient in the HFGMC micromechanical analysis is that all these 

conditions are imposed in the average (integral) sense. 

As a result of the imposition of these conditions, a linear system of algebraic equations is obtained, 

which can be represented in the following form: 

∆=∆ FVK (3)

where the matrix K contains information on the global geometry and instantaneous properties of the 

materials (or composites) within the individual subcells (ΑΒΓ), and the displacement vector increment,

∆V, contains the unknown displacement coefficients )(

)(

ΑΒΓ ∆Wi lmn , which appear on the right-hand side of 

Equation (2). The vector ΔF contains information on the applied average strain increments ε∆ ij . The 

solution of Equation (3) enables the establishment of the following localization relation which expresses 

the average strain increments ΑΒΓ)( ε∆ ij in the global subcell (ΑΒΓ) to the externally applied average 

strain increments ε∆ ij in the form, 

kl ij ijkl A ε∆=ε∆ ΑΒΓ ΑΒΓ)()( (4)

where ΑΒΓ)(

ijkl A are the instantaneous global strain concentration tensor components, of the subcell (ΑΒΓ). 

The final form of the effective incremental constitutive law of the multiphase composite, which 

relates the average stress increments σ∆ ij and strain increments ε∆ kl , is established as follows: 
NASA/TM—2015-218838 6 

kl ij ijkl C∗ ε∆=σ∆ (5)

In this equation ∗

ijkl C are components of the instantaneous effective global stiffness tensor, which are 

given by,

∑ ∑ ∑

Α Β Γ

=Α =Β =Γ

ΑΒΓ ΑΒΓ

ΓΒΑ ∗ =

N N N

ijkl ijpq pqkl ACLHD

DHL

C

1 1 1

1 )()( (6)

Next, the components of the global instantaneous stress concentration tensor, ΑΒΓ)(

ijkl B , which relate the 

average stress increments in the global subcell, ΑΒΓ)( σ∆ ij , to the average (global) stress increments, σ∆ ij , 

are determined. By combining Equations (1) and (4), the subcell stresses are given by, 

ij =σ∆ ijpq AC pqkl ε∆ kl

ΑΒΓ ΑΒΓ ΑΒΓ)()()( (7)

Then, using Equation (5), one obtains, 

ij ijkl kl B σ∆=σ∆ ΑΒΓ ΑΒΓ)()( (8)

where, 

ΑΒΓ ΑΒΓ ΑΒΓ ∗ = ijkl ijpq pqrs rskl SACB

)()()( (9)

and ∗

rskl S are components of the effective instantaneous global compliance tensor, which are the 

components of the inverse of the instantaneous global effective stiffness tensor, C*

, see Equation (6). 

Local Scale Analysis 

As stated above and shown in Figure 2, the global subcells can be occupied by periodic composite 

materials, whose responses are determined by the micromechanical analysis of local RUCs. Therefore the 

global scale analysis subcell quantities are the RUC quantities in the local scale analysis. As such, the 

global subcell incremental constitutive equation (1) serves as the RUC constitutive equation for the local 

scale. The local subcell incremental nonlinear elastic constitutive equation is given by, 

ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ αβγ),(),(),( σ∆ = ε∆ ij ijkl kl C (10)

where the indices (ΑΒΓ) indicate the global subcell and the indices (αβγ) indicate the local subcell. As in 

the global scale analysis, the increment of the displacement vector ΑΒΓ αβγ ),( ∆ i u in each local subcell is 

expanded into quadratic forms in terms the coordinates ),,(

),(

3

),(

2

),(

1

ΑΒΓ α ΑΒΓ β ΑΒΓ γ yyy centered in the 

local subcell, as follows, 
NASA/TM—2015-218838 7 

),(

)002(

2

2),( ),(

3

),(

)020(

2

2),( ),(

2

),(

)200(

2

2),( ),(

1

),(

)001(

),(

3

),(

)010(

),(

2

),(

)100(

),(

1

),(

)000(

)()(),(

4

3

2

1

4

3

2

1

4

3

2

1

ΑΒΓ β ΑΒΓ β ΑΒΓ αβγ ΑΒΓ γ ΑΒΓ γ ΑΒΓ αβγ

ΑΒΓ γ ΑΒΓ αβγ ΑΒΓ α ΑΒΓ α ΑΒΓ αβγ

ΑΒΓ αβγ ΑΒΓ ΑΒΓ αβγ ΑΒΓ α ΑΒΓ αβγ ΑΒΓ β ΑΒΓ αβγ

∆

















∆ + − 















+ −

∆

















∆+ + −

∆ ε∆= ∆+ ∆+ ∆+

i i

i i

i i i

p

i ij j

w

l

w y

h

y

w

d

wy y

u wY wy wy

(11)

where p = Α,Β,Γ for j = 1,2,3, respectively, and d(ΑΒΓ,α), h(ΑΒΓ,β), and l(ΑΒΓ,γ)are the dimensions of the local 

subcell (αβγ) within the global subcell (ΑΒΓ). The unknown local terms ),(

)(

ΑΒΓ αβγ ∆ i lmn w are determined 

from the fulfillment of the equilibrium conditions, the periodic boundary conditions, and the interfacial 

continuity conditions of displacements and tractions between local subcells. Again, all of these conditions 

are imposed in the average sense. 

The application of these conditions lead to the local system of linear algebraic equations, which can 

be written in the following form: 

∆=∆ fvk (12)

where the matrix k contains information on the local geometry and instantaneous properties of the 

materials within the individual subcells (αβγ), and the local displacement vector increment, ∆v, contains 

the unknown local displacement coefficients ),(

)(

ΑΒΓ αβγ ∆wi lmn , which appear on the right-hand side of 

Equation (11). The vector ∆f contains information on the global subcell average strain increments 

ΑΒΓ)( ε∆ ij . The solution of Equation (12) enables the establishment of the following localization relation 

which expresses the average local strain increments ΑΒΓ αβγ ),( ε∆ ij in the local subcell (αβγ) to the global 

subcell average strain increments 

ΑΒΓ)( ε∆ kl in the global subcell (ΑΒΓ),

ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ)(),(),( ε∆ = ε∆ ij ijkl kl A (13)

where ΑΒΓ αβγ ),(

ijkl A are the instantaneous strain concentration tensor components of the local subcell (αβγ) 

within the global subcell (ΑΒΓ). With the established ΑΒΓ αβγ ),(

ijkl A , the local RUC effective instantaneous 

stiffness tensor, which corresponds to the global subcell effective instantaneous stiffness tensor, 

ΑΒΓ)(

ijkl C , 

in Equation (1), can be written as, 

∑ ∑ ∑

α β γ

=α =β =γ

ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ α ΑΒΓ β ΑΒΓ γΑΒΓ =

N N N

ijkl ijpq pqkl Clhd A

dhl

C

1 1 1

),(),( ),(),(),( )( 1

(14)

As before, the components of the local instantaneous stress concentration tensor, ΑΒΓ αβγ ),(

ijkl B , which relate 

the average stress increments in the local subcell, ΑΒΓ αβγ ),( σ∆ ij , to the global subcell average stress
NASA/TM—2015-218838 8 

increments, 

ΑΒΓ)( σ∆ ij , are determined. By combining Equations (10) and (13), the local subcell stresses 

are given by,

ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ)(),(),(),( σ∆ = ε∆ ij ijpq pqkl kl C A (15)

Then, using Equation (1), one obtains, 

(ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ)(),(), σ∆ = σ∆ ij ijkl kl B (16)

where, 

(ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ)(),(),(), = ijkl ijpq pqrs rskl B C SA (17)

and ΑΒΓ)(

rskl S are components of the effective instantaneous global subcell compliance tensor, which 

are the components of the inverse of the instantaneous global subcell stiffness tensor, C(ΑΒΓ)

, see 

Equation (14). 

Therefore, starting at the local scale, given the properties of the local-scale constituents (e.g., fiber 

and matrix) and their arrangement in the RUC, the local HFGMC provides the incremental system of 

algebraic equations (12), from which the instantaneous strain concentration tensor, A(ΑΒΓ,αβγ)

, is 

determined. Equation (14) then provides the instantaneous effective stiffness tensor, C(ΑΒΓ)

, of each

composite material on the local scale . At the global scale, these materials are located within the global 

subcells (ΑΒΓ). Given the arrangement of these composite materials within the global RUC, the global 

HFGMC analysis provides the incremental system of algebraic equations (3), from which the global 

instantaneous strain concentration tensor, A(ΑΒΓ)

, is determined. Then, the effective instantaneous stiffness 

tensor of the entire multiscale material, C*

, can be determined from Equation (6). Given external applied 

loading at the global scale, the stress and strain fields throughout the multiscale material (at both scales) 

can be obtained through the localization equations, Equations (4) and (8) for the global scale and 

Equations (13) and (16) for the local scale.

Initial Yield and Damage Envelopes—Methodology

In the present section, the methodology used to generate initial yield and damage envelopes within 

the multiscale HFGMC analysis is described. The yield surfaces are based on the traditional von Mises 

criteria, as defined by the equivalent stress in the local subcell, given by, 

Y eq ij ij σ σσ= = ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ αβγ ),(),(),( ˆˆ 2

3 (18)

where ˆ 3 ),(),(),( σij σ= ij kk δσ− ij

ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ αβγ are the subcell deviatoric stresses at the local scale, δij

is the Kronecker delta, and Y is the yield stress of the material in the local subcell in simple tension. 

The damage initiation criterion, suggested by Lemaitre and Chaboche (1990), is used to express the 

critical strain energy release rate associated with loss of stiffness in brittle materials, which are typically 

highly dependent on the state of triaxial stress, σh = σkk/3. The average triaxiality function (Lemaitre, 

2001) in the local subcell is given by, 
NASA/TM—2015-218838 9 

( ) ( )

2

),(

),(

),( ),( ),( 1 213

3

2

















σ

σ

ν+= ν−+ ΑΒΓ αβγ

ΑΒΓ αβγ

ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ αβγ

eq

h Rv (19)

where ν(ΑΒΓ,αβγ) is the Poisson’s ratio of the isotropic local subcell and ΑΒΓ αβγ ),( σh is the local subcell 

average hydrostatic pressure. The damage initiation criterion is then given by, 

σeq Rv σ= cr ΑΒΓ αβγ ΑΒΓ αβγ ),(),( (20)

Where σcr is a material parameter that specifies the stress at which damage initiates for a uniaxial tensile

test. 

To generate a yield or damage initiation envelope in a given global scale stress plane, for example, 

σ11 – σ22 , the global stress component increments are expressed as,

sin α=σ∆ cosα=σ∆ 11 22 R R (21)

where R is the radial distance from the origin to a point located on the initial yield or damage surface in the 

σ11 – σ22 stress plane and α is the corresponding polar angle. For a given α, the global subcell stress 

increments, ΑΒΓ)( σ∆ ij , are readily determined (in terms of R) by substituting Equation (21) into Equation (8), 

=σ∆ +α α

ΑΒΓ ΑΒΓ ΑΒΓ sin cos )(

22

)(

11

)( RBRB ij ij ij (22)

Substituting this equation into the local stress concentration equation (16), 

σ∆ = ( +α α) ΑΒΓ αβγ ΑΒΓ αβγ ΑΒΓ ΑΒΓ sin cos )(

22

)(

11

),(),(

ij ijkl kl kl BR B B (23)

Due to initial linearity, the yield or damage initiation envelope can be reached in just one increment, 

therefore, the local subcell stresses can be considered as equivalent to their increments. Then the local 

subcell stresses given by Equation (23), which are proportional to R, can be utilized in either initiation 

criterion, Equation (18) or (20), thus providing the value of R/Y or R/σcr corresponding to yield or damage 

initiation at the specified polar angle α. 

Modeling of Fiber Clustering

Herein, the fiber clustering is modeled using two zones with two different volume fractions, which 

may be arranged arbitrarily. To this end, let the fiber reinforced composite material be composed of two 

distinct regions, zone A and zone B, which are distributed within the composite in an arbitrary manner, 

see Figure 1. Note that the number of zones is arbitrary, but herein we have chosen to use two to 

maximize the effect of clustering. Following Shi et al. (2004), the volume fractions of the fibers within 

each zone are denoted as A

f c and B

f c . Let V be the total volume of the composite and VA

 the volume of 

zone A. In addition, let Vf and A

f V be the total volume of fibers in the composite and within zone A, 

respectively. Consequently, the overall (average) volume fraction of the fibers in the composite is,
NASA/TM—2015-218838 10

V

V

c f

f = (24)

Further, the volume fraction of zone A within the composite is given by, 

V

V

c

A

A = (25)

and the fraction of the fibers that are located in zone A is given by, 

f

A

f

V

V

=ζ (26)

These equations imply that within zone A, the fiber volume fraction is, 

ζ=

A

A f

f c

c

c (27)

and within zone B, the fiber volume fraction is, 

( )

A

f B

f c

c

c −

ζ− = 11

(28)

Note that, when cA = ζ, the fibers are uniformly distributed in the composite, and thus no clustering is 

present. When cA = 1, the entire composite is occupied by zone A, whereas, when cA = 0, the entire 

composite is occupied by zone B, with both of these conditions also representing no clustering. When 

ζ = 1, all fibers are located in zone A, while zone B is pure matrix, which represents complete clustering.

Likewise, when ζ = 0, all fibers are located in zone B, while zone A is pure matrix, which also represents 

complete clustering. Therefore, for a given segregation of the composite into zones A and B, the degree of 

fiber clustering can be fully described by the values of cA and ζ, both of which vary between 0 and 1. 

For a given composite overall volume fraction, cf, and a given volume fraction of zone A within the 

composite, cA, the above equations impose limits on the variation ζ. In order for the fiber volume fractions 

of zones A and B, A

f c and Bf c , to remain less than 1, the maximum and minimum values of ζ are given by,

f

A

f

Af

c

c

c

cc

=ζ

−+

min =ζ max ,

1

(29)

Results and Discussion

The effect of fiber clustering could be investigated by arbitrarily selecting the global scale RUC and 

choosing various values of cA and ζ. Here, on the other hand, a micrograph of the 0° ply within a continuous 

fiber carbon/epoxy composite, see Figure 3, has been used to provide a more realistic global scale RUC. 

The overall fiber volume fraction in this micrograph, determined via image pixel analysis is cf = 0.457. To 

assist in dividing the composite into zones A and B, a 10 by 16 square grid was placed over the micrograph, 
NASA/TM—2015-218838 11

and the fiber volume fraction of each square was determined via image pixel analysis, see Table 1. Then, the 

grid squares could be distributed to zones A and B based on a cutoff value of the fiber volume fraction. That 

is, given a cutoff value for the fiber volume fraction (e.g., 0.55), all grid squares above this value were 

assigned to zone A, while the remaining squares were assigned to zone B. While the choice of cutoff value 

is arbitrary, herein three such values (0.55, 0.45, 0.38) were chosen to investigate the effect of fiber 

clustering, see Figure 4. Note that, once a cutoff value is chosen, the value of cA is determined as the fraction 

of shaded (zone A) squares in Figure 4. The corresponding values of cA are 0.26, 0.5, and 0.74, respectively.

The value of ζ, on the other hand, must be specified as constrained by Equation (29). Thus, for cA = 0.26, 

0 ≤ ζ ≤ 0.56, for cA = 0.5, 0 ≤ ζ ≤ 1, and for cA = 0.74, 0.44 ≤ ζ ≤ 1. Once ζ was chosen, the fiber volume 

fractions within zones A and B were calculated from Equations (27) and (28). 

Figure 3.—Micrograph from the 0° ply in an IM7/MTM-45 

carbon/epoxy composite. The overall fiber volume fraction, 

based on image pixel analysis, is 0.457.

TABLE 1.—FIBER VOLUME FRACTIONS OF THE 10 BY 16 SQUARE GRID THAT DIVIDES THE MICROGRAPH OF 

FIGURE 3. THE SHADED ENTRIES ARE THOSE WITH FIBER VOLUME FRACTIONS ABOVE 0.55, WHICH 

CORRESPOND TO THOSE SHADED IN FIGURE 4(A)

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

1 0.648 0.507 0.715 0.701 0.609 0.466 0.562 0.485 0.573 0.449 0.411 0.431 0.328 0.577 0.271 0.267

2 0.509 0.512 0.435 0.506 0.577 0.575 0.706 0.604 0.662 0.714 0.592 0.645 0.638 0.650 0.382 0.417

3 0.447 0.473 0.357 0.465 0.518 0.483 0.559 0.469 0.548 0.537 0.447 0.387 0.428 0.482 0.450 0.446

4 0.380 0.490 0.554 0.252 0.487 0.381 0.537 0.611 0.586 0.532 0.404 0.362 0.516 0.314 0.373 0.454

5 0.375 0.515 0.612 0.627 0.674 0.602 0.559 0.573 0.573 0.573 0.652 0.468 0.486 0.305 0.362 0.271

6 0.278 0.384 0.366 0.404 0.515 0.692 0.699 0.652 0.654 0.566 0.570 0.642 0.424 0.310 0.234 0.299

7 0.188 0.182 0.253 0.110 0.419 0.438 0.347 0.359 0.421 0.598 0.489 0.409 0.306 0.342 0.200 0.317

8 0.068 0.341 0.507 0.373 0.294 0.507 0.395 0.242 0.359 0.341 0.385 0.386 0.436 0.369 0.464 0.506

9 0.298 0.414 0.382 0.386 0.385 0.465 0.596 0.520 0.455 0.565 0.462 0.374 0.361 0.440 0.406 0.385

10 0.446 0.482 0.389 0.410 0.418 0.442 0.416 0.472 0.451 0.439 0.330 0.500 0.330 0.534 0.597 0.462
NASA/TM—2015-218838 12

Figure 4.—Division of the micrograph (Figure 3) into zones A and B based on a cutoff value of the fiber volume 

fraction. (a) Cutoff fiber volume fraction = 0.55 and cA = 0.26, (b) Cutoff fiber volume fraction = 0.45 and cA = 0.5

(c) Cutoff fiber volume fraction = 0.38 and cA = 0.74.

The microstructures shown in Figure 4 were employed as the global RUCs in the multiscale HFGMC 

analyses. At the local scale, 3 by 3 RUCs were employed to maximize efficiency, where the fiber volume 

fraction in zones A and B were adjusted to obtain the desired local fiber volume fractions, see Figure 5, 

Figure 6, and Figure 7. Note that square fiber packing was assumed at the local scale throughout this 

investigation. Table 2 provides a summary of the analysis cases considered associated with the 

micrograph in Figure 3. 

The transversely isotropic carbon fiber properties used within the local scale RUCs are as 

follows: axial Young’s modulus = 286 GPa, transverse Young's modulus = 12.4 GPa, axial shear 

modulus = 20 GPa, axial Poisson’s ratio = 0.29, and transverse Poisson's ratio = 0.29. The isotropic epoxy 

matrix properties considered are: Young’s modulus = 5.2 GPa and Poisson's ratio = 0.35. Throughout this 

study, the continuous fibers are oriented in the one-direction. 

(a) (b)

(c)
NASA/TM—2015-218838 13

Figure 5.—Multiscale material distribution for cA = 0.26 with (a) ζ = 0.5 and (b) ζ = 0.1. 

Figure 6.—Multiscale material distribution for cA = 0.5 with (a) ζ = 0.74 and (b) ζ = 0.24. 

(a) (b)

(a) (b)
NASA/TM—2015-218838 14

Figure 7.—Multiscale material distribution for cA = 0.74 with (a) ζ = 0.89 and (b) ζ = 0.48. 

TABLE 2.—SUMMARY OF THE FIBER CLUSTERING PARAMETER CASES CONSIDERED

FOR THE COMPOSITE MICROGRAPH SHOWN IN FIGURE 3 

Cut off fiber volume 

fraction

cA ζ A

f c Bf c

Used to divide squares in 

micrograph into 

zones A and B

Volume fraction of 

zone A within the 

composite

Fraction of the fibers 

that are located in 

zone A

Fiber volume fraction 

within zone A

Fiber volume fraction 

within zone B

0.55 0.26 0.5 0.89 0.31

0.1 0.18 0.55

0.45 0.5 0.74 0.68 0.24

0.24 0.22 0.7

0.38 0.74 0.89 0.55 0.2

0.48 0.3 0.93

Effective Elastic Properties

The effect of fiber clustering on the effective properties can be extracted from the effective stiffness 

tensor, C*

. In fact, the determination of C* does not require the coupled multiscale approach employed 

herein. Rather, by homogenizing zones A and B independently, and then employing the resulting stiffness 

tensors in a separate, uncoupled, homogenization procedure, C*

 could also be determined. Here, however, 

the fully coupled approach was used. It is emphasized that the initial yield and damage envelopes 

presented in the next section, based on constituent properties, necessitate the use of a coupled multiscale 

approach as the fields at the local scale are required.

Figure 5 shows the microstructural details of the case where cA = 0.26. This figure illustrates the 

global and local RUC geometries used to consider the cases of ζ = 0.5 and ζ = 0.1. For ζ = 0.5, one half 

of the fibers are located in zone A, which is considerably smaller than zone B. Therefore, the volume 

fraction in zone A is quite high ( A

f c = 0.89), whereas Bf c = 0.31 (maintaining the overall cf = 0.457). In 

(a) (b)
NASA/TM—2015-218838 15

contrast, for ζ = 0.1, only 10 percent of the fibers are located in zone A. Thus the fiber volume fractions 

of zone A and zone B are reversed, with A

f c = 0.18, whereas Bf c = 0.55. Both cases represent deviations 

from uniform distribution of fibers in the composite; i.e., clustering. The additional cases where cA = 0.5 

and 0.74 are depicted in Figure 6 and Figure 7. 

Figure 8 shows the normalized effective axial and transverse (in both directions) Young's moduli (E11, 

E22, and E33) and the normalized effective axial shear modulus (G12) of the composite for different values 

of cA and ζ. For each of the three cases in Figure 4, two values of the parameter ζ, which indicates the 

fraction of fibers located in zone A, are considered: one which provides higher fiber volume fraction in 

zone A and the other that provides higher fiber volume fraction in zone B. The results are normalized by 

the uniform fiber distribution case, which can be obtained by setting cA = ζ. Clearly, the impact of fiber 

clustering on the effective properties of the composite is minor. The maximum variation, which occurred 

for G12, is approximately 6 percent versus the uniform fiber distribution case. This validates the traditional 

approach of predicting composite effective elastic properties by considering a uniform fiber distribution.

It is interesting to note that, while clustering resulted in slight reductions in the effective transverse 

Young's moduli, E22 and E33, it resulted in a slight increase in the effective axial shear modulus, G12. Of 

course, the comparison is made to the case of uniform fiber distribution, which herein assumes square 

fiber packing.

Initial Yield and Damage Envelopes 

In this section, the initial yield and damage envelopes are presented for several global scale stress 

planes while varying the fiber clustering. A parametric study is performed by varying the value of cA as 

shown in Figure 4. As in the case of the effective elastic properties above, for each of the three cases in 

Figure 4, two values of the parameter ζ are considered: one which provides higher fiber volume fraction 

in zone A and the other that provides higher fiber volume fraction in zone B. Details are shown in 

Figure 5 to Figure 7. Also shown in the initial envelopes for comparison is corresponding case with 

uniform fiber distribution. 

Figure 8.—The effect of fiber clustering on the normalized effective elastic properties of the composite.
NASA/TM—2015-218838 16

Figure 9 shows the initial yield and damage envelopes for cA = 0.26 for uniform fiber distribution 

and clustering with ζ = 0.5 and ζ = 0.1. Note that the stress magnitudes are normalized with respect to 

the yield stress (Ym) or critical damage stress (σcr-m) of the isotropic matrix in simple tension, see 

Equations (18) and (20). Further, the predictions of the initial envelopes are based on the yielding and 

damage of the first matrix subcell within the local scale RUCs. That is, the fiber is not permitted to yield 

or damage.

Comparing the initial yield envelopes with their corresponding initial damage envelopes, the effect of 

the inclusion of hydrostatic stress in the damage criterion is seen as the damage envelopes are in all cases 

smaller than the corresponding yield envelopes. This is particularly true in the σ22-σ33 envelopes, wherein 

the global hydrostatic stress is maximized in quadrants I and III. Clearly, the effect of clustering is 

significant since the envelopes that include this effect are always inside the envelopes associated with 

uniform fiber distribution. It should be emphasized that these initiation envelopes may or may not be 

representative of composite ultimate strength. In brittle systems (e.g., ceramic matrix composites), 

however, they would be expected to be more representative of final failure. In addition, the ζ = 0.5 

envelopes are always inside the ζ = 0.1 envelopes. This is likely caused by the very high fiber volume 

fraction of zone A ( A

f c = 0.89) in the case of ζ = 0.5, which results in high local scale stress 

concentrations in the matrix. 

It is interesting that initiation of yielding and damage may occur in either zone A or zone B at each 

point on each envelope, thus diamond symbols are included in Figure 9 to indicate the portions of each 

envelope where initiation occurred in zone A. Transitions between zone A initiation and zone B initiation 

occur in many of the presented envelopes, while some exhibit initiation only in zone B (solid lines with 

no symbols). Obviously, for the case of uniform fiber distribution, there is no distinction between zones A 

and B, thus such transitions are not applicable. 

Examining the initial σ11-σ22 yield and damage envelopes reveals that, as expected, the effect of the 

continuous fibers is dominant as there is approximately a 25:1 difference in initiation stress for pure σ11 

versus pure σ22. Also, the effect of clustering is negligible under pure axial loading. Also of note is the 

fact that the damage initiation always occurs in zone B irrespective of the value of ζ. In contrast, this 

parameter has a pronounced effect on the initiation location in the case of yielding. 

Considering biaxial loading in the transverse directions (σ33-σ22), the effect of the fiber is minimized.

Here, all of the envelopes are nearly symmetric, and the locations of the transitions between zone A and 

zone B initiation are qualitatively similar between the yield and damage surfaces. 

The effect of clustering is most evident in the initial σ12-σ22 yield and damage envelopes. The predicted 

envelopes appear to be concentric ellipses whose semimajor axes are dictated by the clustering. For 

ζ = 0.5, the yield and damage envelopes exhibit similar transitions from zone A to zone B initiation. 

However, for ζ = 0.1, such transitions are present only in the yield envelope. It is interesting to note that 

while the initial yield and damage points under pure axial shear (which are based on the local fields) are 

reduced by the presence of clustering, the effective axial shear moduli shown in Figure 8 (which are based 

on the two-scale homogenized stiffness tensor) show an increase due to clustering.

The microstructural details in the case of cA = 0.5 and ζ = 0.74 and ζ = 0.24 are shown in Figure 6. 

The corresponding initial yield and damage envelopes are shown in Figure 10. Similarly, for the case of 

cA = 0.74 and ζ = 0.89 and ζ = 0.48, the microstructural details are given in Figure 7, and the associated 

initial envelopes are given in Figure 11. A comparison of Figure 9, Figure 10, and Figure 11, in which the 

volume fraction of zone A, cA, increases from 0.26 to 0.5 to 0.74, reveals that, while all representations of 

clustering provide smaller yield envelopes compared to the uniform case, the ordering of high ζ and low ζ

envelopes switches. In Figure 9, the higher ζ envelopes are smallest, in Figure 10 the higher and lower ζ 
NASA/TM—2015-218838 17

Initial Yield Envelopes Initial Damage Envelopes

Figure 9.—Initial yield and damage envelopes for cA = 0.26 associated with Figure 5. The symbols on certain 

envelopes indicate yield/damage initiation within zone A, whereas, when no symbol is present, failure initiated in 

zone B.

-30

-20

-10

0

10

20

30

σ-2 -1 0 1 2 11/Y

m

σ22/Ym

 zeta = 0.5

 uniform

 zeta=0.1

-30

-20

-10

0

10

20

30

-2 -1 0 1 2

σ11

/

σcr-m

σ22/σcr-m

 zeta = 0.5

 uniform

 zeta = 0.1

-4

-2

0

2

4

-4 -2 0 2 4 σ33/Y

m

σ22/Ym

 zeta = 0.5

 uniform

 zeta = 0.1

-4

-3

-2

-1

0

1

2

3

4

-4 -2 0 2 4

σ33

/

σcr-m

σ22/σcr-m

 zeta = 0.5

 uniform

 zeta = 0.1

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5

σ12/Y

m

σ22/Ym

 zeta = 0.5

 uniform

 zeta = 0.1

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5

σ12/σcr-m

σ22/σcr-m

 zeta = 0.5

 uniform

 zeta = 0.1
NASA/TM—2015-218838 18

Initial Yield Envelopes Initial Damage Envelopes

Figure 10.—Initial yield and damage envelopes for cA = 0.5 associated with Figure 6. The symbols on certain 

envelopes indicate yield/damage initiation within zone A, whereas, when no symbol is present, failure initiated in 

zone B.

-30

-20

-10

0

10

20

30

σ

-2 -1 0 1 2

11/Y

m

σ22/Ym

 zeta = 0.74

 uniform

 zeta=0.24

-30

-20

-10

0

10

20

30

-2 -1 0 1 2

σ11

/

σcr-m

σ22/σcr-m

 zeta = 0.24

 uniform

 zeta = 0.74

-4

-2

0

2

4

σ-4 -2 0 2 4 33/Y

m

σ22/Ym

 zeta = 0.74

 uniform

 zeta = 0.24

-4

-3

-2

-1

0

1

2

3

4

-4 -2 0 2 4 σ33

/

σcr-m

σ22/σcr-m

 zeta = 0.74

 uniform

 zeta = 0.24

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5 σ12/Y

m

σ22/Ym

 zeta = 0.74

 uniform

 zeta = 0.24

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5

σ12/σcr-m

σ22/σcr-m

 zeta = 0.74

 uniform

 zeta = 0.24
NASA/TM—2015-218838 19

Initial Yield Envelopes Initial Damage Envelopes

Figure 11.—Initial yield and damage envelopes for cA = 0.74 associated with Figure 7. The symbols on certain 

envelopes indicate yield/damage initiation within zone A, whereas, when no symbol is present, failure initiated in 

zone B. Note that in the σ11-σ22 envelopes for ζ = 0.48, under pure σ11 loading, yield and damage initiated in zone B 

(these two points on these envelopes have no symbol).

-30

-20

-10

0

10

20

30

σ-2 -1 0 1 2 11/Y

m

σ22/Ym

 zeta = 0.89

 uniform

 zeta=0.48

-30

-20

-10

0

10

20

30

-2 -1 0 1 2

σ11

/

σcr-m

σ22/σcr-m

 zeta = 0.89

 uniform

 zeta = 0.48

-4

-2

0

2

4

-4 -2 0 2 4 σ33/Y

m

σ22/Ym

 zeta = 0.89

 uniform

 zeta = 0.48

-4

-2

0

2

4

-4 -2 0 2 4

σ33

/

σcr-m

σ22/σcr-m

 zeta = 0.89

 uniform

 zeta = 0.48

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5 σ12/Y

m

σ22/Ym

 zeta = 0.89

 uniform

 zeta = 0.48

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5

σ12/σcr-m

σ22/σcr-m

 zeta = 0.89

 uniform

 zeta = 0.48
NASA/TM—2015-218838 20

envelopes are very close, and in Figure 11, the lower ζ envelopes are smallest. This observation is 

reasonable since the smallest envelopes are always associated with the higher local fiber volume fraction 

being located in the zone with the lower volume fraction within the global RUC. This situation represents 

the greater degree of fiber clustering and leads to higher stress concentrations in the composite, and thus 

contracted initial envelopes. It should be noted that, as before, the transitions of yield and damage 

initiation from zones A to B are denoted in Figure 10 and Figure 11. 

To examine the effect of cA on the initial yield and damage envelopes in a consistent way, the actual ζ

value associated with each of the discretizations shown in Figure 4 can be calculated. This is done by 

rearranging Equation (27), 

f

A A

act f c

c =ζ c (30)

where A

f c , the average volume fraction in zone A, is calculated for each cA value by simply averaging the 

fiber volume fractions of the grid squares in zone A. For example, for cA = 0.26, A

f c is equal to the 

average of the shaded entries in Table 1. For cA = 0.26, ζact = 0.35 ( A

f c = 0.62 and Bf c = 0.40); for 

cA = 0.5, ζact = 0.61 ( A

f c = 0.56 and Bf c = 0.36); for cA = 0.74, ζact = 0.83 ( Af c = 0.51 and Bf c = 0.30). 

The effect of cA on the predicted composite initial yield and damage envelopes of the carbon/epoxy 

composite, where ζ = ζact, is shown in Figure 12. Clearly, the variation among the envelopes is less than 

that shown in Figure 9 to Figure 11. This is expected as, for all three cases of cA = 0.26, 0.5, and 0.74, the 

ζact value is between the values of ζ used in Figure 9 to Figure 11. However, this also indicates that the 

choice of cA used to discretize the microstructure is not critical.

In order to illustrate the effect of material system on the impact of fiber clustering, the initial yield 

and damage envelopes of an E-glass/epoxy composite have also been examined. The considered isotropic 

elastic properties of the E-glass fibers are: E = 73 GPa and ν =0.22, while the same epoxy matrix is 

considered. The mismatch in transverse Young's modulus between the fiber and matrix is known to have 

a significant effect on the local fields within composites. For the carbon/epoxy composite considered 

above, this mismatch is 12.4/5.2 = 2.4, whereas for the glass/epoxy composite, it is 73/5.2 = 14. Note that 

the identical carbon/epoxy microstructure (Figure 3) and idealizations into zones (Figure 4) were used.

Figure 13 shows the effect of cA on the predicted composite initial yield and damage envelopes of the 

E-glass/epoxy composite, where ζ = ζact. Compared to the carbon/epoxy composite (Figure 12), it is clear 

that the effect of clustering on the initial yield and damage envelopes is much greater in the case of the 

E-glass/epoxy composite. The E-glass/epoxy’s greater mismatch in the fiber-matrix transverse Young’s 

modulus leads to higher stress concentrations due to the presence of clustering, which are manifested in 

the predicted envelopes. Similarly, there are now some noticeable variations in the envelopes based on the 

choice of cA, but, as in the case of the carbon/epoxy composite, this difference is significantly less than 

the difference between the uniform case and any of the clustered cases.
NASA/TM—2015-218838 21

Initial Yield Envelopes Initial Damage Envelopes

Figure 12.—Comparison of initial yield and damage envelopes for the carbon/epoxy composite, for different values 

of cA with ζ = ζact. 

-30

-20

-10

0

10

20

30

-2 -1 0 1 2 σ11/Y

m

σ22/Ym

cA = 0.26

cA = 0.5

cA = 0.74

 uniform

-30

-20

-10

0

10

20

30

-2 -1 0 1 2

σ11

/

σcr-m

σ22/σcr-m

cA = 0.26

cA = 0.5

cA = 0.74

 uniform

-4

-3

-2

-1

0

1

2

3

4

-4 -3 -2 -1 0 1 2 3 4 σ33/Y

m

σ22/Ym

cA = 0.26

cA = 0.5

cA = 0.74

 uniform

-4

-3

-2

-1

0

1

2

3

4

-4 -3 -2 -1 0 1 2 3 4

σ33

/

σcr-m

σ22/σcr-m

cA = 0.26

cA = 0.5

cA = 0.74

 uniform

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5 σ12/Y

m

σ22/Ym

cA = 0.26

cA = 0.5

cA = 0.74

 uniform

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5

σ12

/

σcr-m

σ22/σcr-m

cA = 0.26

cA = 0.5

cA = 0.74

 uniform
NASA/TM—2015-218838 22

Initial Yield Envelopes Initial Damage Envelopes

Figure 13.—Comparison of initial yield and damage envelopes for the E-glass/epoxy composite, for different values 

of cA with ζ = ζact. 

-8

-6

-4

-2

0

2

4

6

8

-2 -1 0 1 2 σ11/Y

m

σ22/Ym

uniform

cA = 0.26

cA = 0.5

cA = 0.74

-8

-6

-4

-2

0

2

4

6

8

-2 -1 0 1 2

σ11

/

σcr-m

σ22/σcr-m

uniform

cA = 0.26

cA = 0.5

cA = 0.74

-2.5

-2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

2.5

-2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 σ33/Y

m

σ22/Ym

cA = 0.26

cA = 0.5

cA = 0.74

uniform

-2.5

-2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

2.5

-2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5

σ33

/

σcr-m

σ22/σcr-m

cA = 0.26

cA = 0.5

cA = 0.74

uniform

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 σ12/Y

m

σ22/Ym

cA = 0.26

cA = 0.5

cA = 0.74

uniform

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5

σ12

/

σcr-m

σ22/σcr-m

cA = 0.26

cA = 0.5

cA = 0.74

uniform
NASA/TM—2015-218838 23

Nonlinear Effective Stress-Strain Curves

To examine the effect of fiber clustering on the nonlinear behavior of the composite, the epoxy matrix 

material response was simulated using the nonlinear elastic Ramberg and Osgood (1943) model. The 

uniaxial form of the constitutive relation in this model is given by, 

n

EE 















σ

σσ

+

σ

=ε

0

0 (31)

where E is the Young’s modulus and σ0 and n are two material parameters that can be determined by 

fitting the material nonlinear stress-strain response. A deviatoric multiaxial generalization is given by

1

0

ˆ 2

1 3 −

















σ

σ

σ+δσ ν−σ

ν+

=ε

n eq

ij ij kk ij ij EE E (32)

where ν is the Poisson’s ratio and σeq is the equivalent stress (see Eq. (18)). 

In the implementation of Equation (32) within the multiscale HFGMC model described above, it is 

necessary to cast this equation in its incremental form,

Sijklij σ∆=ε∆ kl (33)

where Sijkl are the components of the instantaneous compliance tensor, which can be obtained from 

manipulation of the incremental form of Equation (32). Then, inverting S, the local subcell constitutive 

equation (10), including the Ramberg-Osgood nonlinear effects, is obtained. The remainder of the 

incremental multiscale HFGMC theory remains the same as presented above. 

The nonlinearity of the epoxy material has been modeled herein by choosing two sets of parameters 

with different degrees of nonlinearity: σ0 = 158 MPa, n = 4 (from Aboudi, 1991) and σ0 = 80 MPa, n = 15 

(fictional values, resulting in greater nonlinearity). The resulting uniaxial stress-strain curves in simple 

tension are shown in Figure 14. 

The resulting effective transverse stress-strain curves of the composite at the global scale are shown 

in Figure 15 for the three values of cA with varying ζ and for both sets of nonlinear parameters. It can be 

readily observed that the effect of clustering on the C/epoxy composite is greater when more matrix 

nonlinearity is present. It would be expected that this effect would continue to increase with decreasing 

matrix instantaneous stiffness. In addition, by comparing the three plots in Figure 15, it can be seen that 

for cA = 0.26, the curve associated with higher ζ has the greater deviation from the uniform case, whereas 

for cA = 0.74, the curve associated with lower ζ has the greater deviation from the uniform case. For 

intermediate value of cA = 0.5, the curves associates with both values of ζ are very similar. As in the case 

of the initial envelopes, the greatest deviation from the uniform results is always associated with the 

higher local fiber volume fraction being located in the zone with the lower volume fraction within the 

global RUC.

Once again, to illustrate the effect of the contrast between the properties of the fibers and the matrix 

on the impact of fiber clustering, the stress-strain response of the E-glass/epoxy composite has been 

examined. The more nonlinear epoxy Ramberg-Osgood parameters, σ0 = 80 MPa, n = 15, were employed 

in the following simulations.
NASA/TM—2015-218838 24

Figure 14.—Uniaxial stress-strain curves of the nonlinear epoxy in 

simple tension, as represented by the Ramberg-Osgood model with 

E = 5.2 GPa, ν = 0.35 and the two sets of nonlinear parameters 

indicated.

Figure 15.—Effective transverse response of the C/epoxy composite using two sets on nonlinear parameters for 

the matrix. (a) cA = 0.26, (a) cA = 0.5, (a) cA = 0.74. 

0

20

40

60

80

100

120

140

160

180

0 0.01 0.02 0.03 0.04 0.05 0.06

σ (MPa)

ε

s0 = 158 MPa

n = 4

s0 = 80 MPa

n = 15

0

50

100

150

200

250

300

350

0 0.01 0.02 0.03 0.04 0.05 0.06

σ22 (MPa)

ε22

uniform

 cA = 0.26, z = 0.5

 cA = 0.26, z = 0.1

s0 = 158 MPa

n = 4

s0 = 80 MPa

n = 15

0

50

100

150

200

250

300

350

0 0.01 0.02 0.03 0.04 0.05 0.06

σ22 (MPa)

ε22

uniform

 cA = 0.5, z = 0.74

 cA = 0.5, z = 0.24

s0 = 158 MPa

n = 4

s0 = 80 MPa

n = 15

0

50

100

150

200

250

300

350

0 0.01 0.02 0.03 0.04 0.05 0.06

σ22 (MPa)

ε22

uniform

 cA = 0.74, z = 0.89

 cA = 0.74, z = 0.48

s0 = 158 MPa

n = 4

s0 = 80 MPa

n = 15

(a) (b)

(c)
NASA/TM—2015-218838 25

Figure 16.—Effective transverse response of the E-glass/epoxy composite compared with the C/epoxy 

composite. (a) cA = 0.26, (a) cA = 0.5, (a) cA = 0.74. 

The results of the transverse tensile simulations of the E-glass/epoxy composite are compared to the 

C/epoxy composite in Figure 16. Obviously, because the E-glass fiber has a higher transverse stiffness, 

the E-glass/epoxy composite also has a higher effective transverse stiffness. In addition, the higher 

transverse mismatch in fiber-matrix properties results in a much greater influence of fiber clustering on 

the effective transverse response. Otherwise, the trend with respect to the parameter ζ is the same in the 

E-glass/epoxy composite as in the C/epoxy system. 

Alternative Microstructure

To investigate the influence of the assumed global microstructure, an alternative global 

microstructure, taken from a different specimen of the same type of composite, was considered. The 

micrograph is shown in Figure 17. As before, a square grid was placed over the micrograph and the local 

fiber volume fraction within each square was determined via image pixel analysis. Because the aspect 

ratio of the present micrograph shown in Figure 17 is different than the previous micrograph shown in 

Figure 3, a grid of 10 x 24 squares was required. The overall fiber volume fraction, cf, determined via 

image pixel analysis of the alternative micrograph is 0.452. A comparison of Figure 5 and Figure 18 

indicates that the two global RUCs are quite different, with the latter having the majority of zone A 

located in the right half of the RUC.

0

50

100

150

200

250

300

350

400

450

500

0 0.01 0.02 0.03 0.04 0.05 0.06

σ22 (MPa)

ε22

uniform

 cA = 0.26, z = 0.5

 cA = 0.26, z = 0.1

E-glass/epoxy

C/epoxy

0

50

100

150

200

250

300

350

400

450

500

0 0.01 0.02 0.03 0.04 0.05 0.06

σ22 (MPa)

ε22

uniform

 cA = 0.5, z = 0.74

 cA = 0.5, z = 0.24

E-glass/epoxy

C/epoxy

0

50

100

150

200

250

300

350

400

450

500

0 0.01 0.02 0.03 0.04 0.05 0.06

σ22 (MPa)

ε22

uniform

 cA = 0.74, z = 0.89

 cA = 0.74, z = 0.48

E-glass/epoxy

C/epoxy

(a) (b)

(c)
NASA/TM—2015-218838 26

Figure 17.—Alternative micrograph from a 0° ply in an IM7/MTM-45 

carbon/epoxy composite. The overall fiber volume fraction, based on 

image pixel analysis, is 0.452.

Figure 18.—Multiscale material distribution for cA = 0.25 with ζ = 0.5, for the alternative composite 

micrograph shown in Figure 17. 

In order to compare to the previous results, a cA = 0.25 was obtained using a fiber volume fraction cutoff of 

0.55. Note that it was not possible to match the previous cA = 0.26 exactly. The resulting global scale RUC is 

shown in Figure 18, along with the local RUCs and local fiber volume fractions corresponding to ζ = 0.5.

A comparison of the initial yield and damage envelopes of the original microstructure with the present 

alternative microstructure is shown in Figure 19. These comparisons are shown for cA = 0.25 or 0.26 and 

ζ = 0.5. As can be seen, although the micrographs used to generate the global RUCs are quite different, the 

predicted initial envelopes under normal loadings are nearly identical. However, for the axial shear case, the 

envelopes show variation under significant shear loading. Thus because of this variation in the shear 

responses, it is not possible to claim that the clustering is fully characterized by the parameters cA and ζ. 

Rather, the chosen global scale microstructure can also have an effect. It should be noted, however, that, 

since only two particular microstructures were considered, the sample size is in no way statistically 

representative. Similarly, no assessment has been made of the statistical relevance of two specific 

microstructures considered in terms of their size and number of fibers included.
NASA/TM—2015-218838 27

A comparison between the global transverse stress-strain curves obtained from the original 

microstructure and the present alternative microstructure is shown in Figure 20, in which cA = 0.25 or 

0.26 and ζ = 0.5. It can be observed that the two responses are again quite similar.

Initial Yield Envelopes Initial Damage Envelopes

Figure 19.—Comparison of initial yield and damage envelopes for cA = 0.26 and ζ = 0.5 from the original (Figure 5) 

and alternative (Figure 18) microstructures.

-30

-20

-10

0

10

20

30

σ

-2 -1 0 1 2

11/Y

m

σ22/Ym

Original Microstructure

Alternative Microstructure -30

-20

-10

0

10

20

30

-2 -1 0 1 2

σ11

/

σcr-m

σ22/σcr-m

Original microstructure

Alternative microstructure

-4

-2

0

2

4

-4 -2 0 2 4 σ33/Y

m

σ22/Ym

Original microstructure

Alternative microstructure -4

-3

-2

-1

0

1

2

3

4

-4 -2 0 2 4

σ33

/

σcr-m

σ22/σcr-m

Original microstructure

Alternative microstructure

-1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5

σ12/Y

m

σ22/Ym

Original microstructure

Alternative microstructure -1.5

-1

-0.5

0

0.5

1

1.5

-1.5 -1 -0.5 0 0.5 1 1.5

σ12/σcr-m

σ22/σcr-m

Original microstructure

Alternative microstructure
NASA/TM—2015-218838 28

Figure 20.—Comparison of the effective transverse response of the 

C/epoxy composite for the original and alternative microstructures 

with cA = 0.26 and ζ = 0.5. The nonlinear epoxy parameters, 

σ0 = 80 MPa, n = 15, were employed.

Conclusions

A multiscale micromechanical approach has been developed for the prediction of the behavior of 

composites exhibiting fiber clustering. The method is able to predict not only effective composite 

properties, but also initial yield and damage envelopes, as well as the nonlinear stress-strain response of 

the composite, due to the availability of local fields. This is possible because this approach is based on the 

establishment of the concentration tensors on the local and global scales. Fiber clustering is characterized 

by two nondimensional parameters, which measure the extent and density of the fiber clusters, along with 

a global RUC determined from actual composite micrographs. This methodology was applied to 

carbon/epoxy and E-glass/epoxy composites. It has been shown that fiber clustering results in a decrease 

in the sizes of the initial yield and damage surfaces as compared to the uniformly distributed case. The 

effect is greater in the case of E-glass/epoxy compared to the baseline carbon/epoxy due to greater 

mismatch in the fiber-matrix transverse Young’s modulus mismatch in the case of the E-glass/epoxy. In

addition, transverse stress-strain simulations indicated that the effect of fiber clustering increases with 

increasing fiber-matrix property mismatch and with increasing nonlinearity of the matrix.

The present method, which can consider three dimensional composite microstructures, could also be 

applied to the simulation of nanocomposites, whose low volume fractions tend to result in significant 

amounts of clustering. Ceramic matrix composites, which typically have brittle matrices and lower fiber 

volume fractions than polymer matrix composites, would also be expected to exhibit significant effects of 

fiber clustering. Further, since HFGMC has been formulated to include inelasticity, fiber-matrix 

debonding (Aboudi et al., 2013), fiber misalignment (Bednarcyk et al., 2014), and progressive damage 

(Bednarcyk et al., 2010; Pineda et al., 2013), these effects can be incorporated at the local scale in the 

present multiscale implementation to determine the impact of their combination with fiber clustering on 

the composite behavior.

0

50

100

150

200

250

0 0.01 0.02 0.03 0.04 0.05 0.06

σ22 (MPa)

ε22

Original microstructure

New microstructure
NASA/TM—2015-218838 29

References

Abhilash, A.S., Joshi, S.P., Mukherjee, A., and Mishnaevsky Jr., L. (2011) Micromechanics of diffusioninduced damage evolution in reinforced polymers. Comp. Sci. Tech. 71, 333-342. 

Aboudi, J., Arnold, S.M., and Bednarcyk, B.A. (2013) Micromechanics of Composite Materials A 

Generalized Multiscale Approach, Elsevier, New York.

Aboudi, J. (1991) Mechanics of Composite Materials A Unified Micromechanical Approach, Elsevier, 

New York.

Acton, K. and Graham-Brady, L. (2010) Elasto-plastic mesoscale homogenization of composite materials. 

J. Engr. Mech. 136, 613-624. 

Al-Ostaz, A. Diwakar, A., and Alzebdeh, K.I. (2007) Statistical model for characterizing random 

microstructure of inclusion-matrix composites. J. Mater. Sci. 42, 7016-7030. 

Baxter, S.C., Graham-Brady, L.L., Xu, X.F. (2005) Modeling the effects of material nonlinearity using 

moving window micromechanics. Int. J. Nonlinear Mech. 40, 351-359. 

Baxter, S.C. and Graham, L.L. (2000) Characterization of random composite materials using the movingwindow generalized method of cells technique. J. Engr. Mech. 126, 389-397. 

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

to multiaxial progressive damage in the constituents. AIAA J. 48, 1367-1378. 

Bednarcyk, B.A., Aboudi, J., and Arnold, S.M. (2014) The effect of general statistical fiber misalignment 

on damage initiation in composites. Composites: Part B, 66, 97-108. 

Bednarcyk, B.A., Mital, S.K., Pineda, E.J., and Arnold, S.M. (2015) Multiscale Modeling of Ceramic 

Matrix Composites, 56th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and 

Materials Conference, January 5-9, Kissimmee, FL.

Ghosh, S., Nowak, Z., and Lee, K. (1997) Quantitative characterization and modeling of composite 

microstructures by Voronoi cells. Acta Mater. 45, 2215-2234. 

Graham-Brady, L.L., Siragy, E.F., Baxter, S.C. (2003) Analysis of heterogeneous composite materials 

using moving-window techniques. J. Engr. Mech. 129, 1054-1064. 

Guild, F.J. and Summerscales, J. (1993) Microstructural image analysis applied to fibre composite 

materials: a review. Composites 24, 383-393. 

Lemaitre, J. (2001) Handbook of Materials Behavior Models, Elsevier, New York. 

Lemaitre, J. and Chaboche, J.L. (1990) Mechanics of Solid Materials, Cambridge University Press, 

Cambridge. 

Liu, K.C. and Arnold, S.M. (2013) Influence of scale specific features on the progressive damage of 

woven ceramic matrix composites (CMCs). CMC 35, 35-65. 

Liu, K.C., Chattopadhyay, A., Bednarcyk, B.A., Arnold, S.M. (2011) Efficient multiscale modeling 

framework for triaxially braided composites using generalized method of cells. Journal of Aerospace 

Engineering 42, 162-170. 

Melro, A.R., Camanho, P.P., and Pinho, S.T. (2008) Generation of random distribution of fibres in longfibre reinforced composites. Comp. Sci. Tech. 68, 2092-2101. 

Pindera, M.-J. and Aboudi, J. (1988) Micromechanical analysis of yielding of metal matrix composites. 

Int. J. Plasticity 4, 195-214. 

Pineda, E.J., Bednarcyk, B.A., Waas, A.M., and Arnold, S.M. (2013) Progressive failure of a 

unidirectional fiber-reinforced composite using the method of cells: discretization objective 

computational results. Int. J. Solids Struct. 50, 1203-1216. 

Pyrz, R. (1994) Quantitative description of the microstructure of composites. Part I: Morphology of 

unidirectional composite systems. Comp. Sci. Tech. 50, 197-208. 

Ramberg, W. and Osgood, W. R. (1943) Description of stress-strain curves by three parameters. 

Technical Note No. 902, National Advisory Committee For Aeronautics, Washington DC. 
NASA/TM—2015-218838 30

Romanov, V., Lomov, S.V., Swolfs, Y., Orlova, S., Gorbatikh, L., and Verpoest, I. (2013) Statistical 

analysis of real and simulated fibre arrangements in unidirectional composites. Comp. Sci. Tech. 87, 

126-134. 

Scalon, J.D., Fieller, N.R.J., Stillman, E.C., and Atkinson, H.V. (2003) Spatial pattern analysis of secondphase particles in composite materials. Mater. Sci. Engr. A356, 245-257. 

Shaffer, M.S.P. and Windle, A.H. (1999) Fabrication and characterization of carbon nanotube/poly(vinyl 

alcohol) composites. Adv. Mater. 11, 937-941. 

Shi, D.-L., Feng, X.-Q., Huang, Y.Y., Hwang, K.-C., and Gao, H. (2004) The effect of nanotube waviness 

and agglomeration on the elastic properties of carbon nanotube-reinforced composites. J. Eng. Mater. 

Tech. 126, 250-257. 

Sorensen, B.F. and Talreja, R. (1993) Effects of nonuniformity of fiber distribution on the thermallyinduced residual stresses and cracking in ceramic matrix composites. Mech. Mater. 16, 351-363. 

Vaughan, T.J. and McCarthy, C.T. (2010) A combined experimental-numerical approach for generating 

statistically equivalent fiber distributions for high strength laminated composite materials. Comp. Sci. 

Tech. 70, 291-297. 

Vigolo, B., Penicaud, A.P., Couloun, C., Sauder, S., Pailler, R., Journet, C., Bernier, P., and Poulin, P. 

(2000) Macroscopic fibers and ribbons of oriented carbon nanotubes. Science 290, 1331-1334. 

Wilding, S.E. and Fullwood, D.T. (2011) Clustering metrics for two-phase composites. Comput. Mater. 

Sci. 50, 2262-2272. 

Wongsto, A. and Li, S. (2005) Micromechanical FE analysis of UD fibre-reinforced composites with 

fibres distributed at random over the transverse cross-section. Composites A 36, 1246-1266. 

Yang, S., Gokhale, A.M., and Shan, Z. (2000) Utility of microstructure modeling for simulation of micromechanical response of composites containing non-uniformly distributed fibers. Acta mater. 48, 

2307-2322. 

Yang, S. Tewari, A., and Gokhale, A.M. (1997) Modeling of non-uniform spatial arrangement of fibers in 

ceramic matrix composites. Acta mater. 45, 3059-3069. 

Zangenberg, J., Larsen, J.B., Østergaard, R.C., and Brøndsted, P. (2012) Methodology for 

characterisation of glass fibre composite architecture. Plastics, Rubber and Composites 41, 187-193. 

Zeman, J. and Sejnoha, M. (2001) Numerical evaluation of effective elastic properties of graphite fiber 

tow impregnated by polymer matrix. J. Mech. Phys. Sol. 49, 69-90.