Evan J. Pineda, Brett A. Bednarcyk, and Steven M. Arnold

Glenn Research Center, Cleveland, Ohio

Effects of Subscale Size and Shape on Global 

Energy Dissipation in a Multiscale Model of a 

Fiber-Reinforced Composite Exhibiting Post-Peak 

Strain Softening Using Abaqus and FEAMAC

NASA/TM—2012-217650

August 2012
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Officer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of significant 

scientific and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question 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, Brett A. Bednarcyk, and Steven M. Arnold

Glenn Research Center, Cleveland, Ohio

Effects of Subscale Size and Shape on Global 

Energy Dissipation in a Multiscale Model of a 

Fiber-Reinforced Composite Exhibiting Post-Peak 

Strain Softening Using Abaqus and FEAMAC

NASA/TM—2012-217650

August 2012

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

Simulia Community Conference (SCC)

sponsored by the Dassault Systèmes

Providence, Rhode Island, May 15–17, 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 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/TM—2012-217650 1

Effects of Subscale Size and Shape on Global Energy Dissipation in 

a Multiscale Model of a Fiber-Reinforced Composite Exhibiting 

Post-Peak Strain Softening Using Abaqus and FEAMAC

Evan J. Pineda, Brett A. Bednarcyk, and Steven M. Arnold

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Abstract

A mesh objective crack band model is implemented in the generalized method of cells (GMC) 

micromechanics model to predict failure of a composite repeating unit cell (RUC). The micromechanics 

calculations are achieved using the MAC/GMC core engine within the ImMAC suite of micromechanics 

codes, developed at the NASA Glenn Research Center. The microscale RUC is linked to a macroscale 

Abaqus/Standard finite element model using the FEAMAC multiscale framework (included in the 

ImMAC suite). The effects of the relationship between the characteristic length of the finite element and 

the size of the microscale RUC on the total energy dissipation of the multiscale model are investigated. A 

simple 2-D composite square subjected to uniaxial tension is used to demonstrate the effects of scaling the 

dimensions of the RUC such that the length of the sides of the RUC are equal to the characteristic length 

of the finite element. These results are compared to simulations where the size of the RUC is fixed, 

independent of the element size. Simulations are carried out for a variety of mesh densities and element 

shapes, including square and triangular. Results indicate that a consistent size and shape must be used to 

yield preserve energy dissipation across the scales.

1.0 Introduction

Often, details at the constituent level drive the failure of a composite material. An effective way to 

incorporate these details into a numerical model is through multiscale analysis. However, careful attention 

must be paid to how failure is modeled at the different scales. If the stress-strain response of the material 

enters the post-peak softening regime and the tangent stiffness tensor looses positive definiteness, the 

results become sensitive to changes in the mesh. Even if mesh objective constitutive models are in place 

at the microscales or macroscales, mesh objectivity is not necessarily maintained across the scales, and 

the resulting multiscale model will still exhibit pathological dependence on the size of the discretizations 

in one or more scales due to the lack of a clearly defined characteristic length(s) that bridges the scales 

(Bažant, 2007).

A multiscale model is presented here which couples the generalized method of cells (GMC) 

micromechanics model (Paley and Aboudi, 1992), to a macroscopic finite element method (FEM) model.

The commercial FEM software Abaqus is used as the macroscale platform (Simulia Corp., 2008).

Whereas, the openly distributed Integrated Multiscale Micromechanics Analysis Code (ImMAC) suite of 

micromechanics codes, developed by the NASA Glenn Research Center (GRC), is utilized to perform the 

microscale calculations (Bednarcyk and Arnold, 2002a; Bednarcyk and Arnold, 2002b). The 

Micromechanics Analysis Code/Generalized Method of Cells (MAC/GMC) core engine, within the 

ImMAC suite, is used to implement and solve the subscale GMC problem, and the FEAMAC multiscale 

framework (Bednarcyk and Arnold, 2006) facilitates the communication between Abaqus and 

MAC/GMC via Abaqus user subroutines.

A mesh objective crack band model (Bažant and Oh, 1983) was implemented at the subcell level 

within the GMC micromechanics model. With the crack band model, the negative, post peak, tangent 

slope of the stress-strain response of the material is scaled so that the strain energy release rate (SERR) of 
NASA/TM—2012-217650 2

the subcell equals the critical fracture toughness of the material when the stress in the subcell has 

diminished to zero. The crack band model at the microscale is used to influence the behavior of the 

macroscale, in the multiscale FEAMAC model. Pineda et al. (2012) previously implemented the crack 

band model within the high fidelity generalized method of cells (HFGMC), and demonstrated the 

objectivity of the HFGMC solution with respect to the subcell meshes used in the repeating unit cell 

(RUC). By extension, it is assumed that mesh objectivity is assured within the GMC framework at the 

microscale, but this does not warrant mesh objectivity at the macroscale.

The objective of this work is to investigate the effect of the size and shape of the subscale RUC, 

related to the macroscale finite element, on the preservation of energy in the multiscale model. Thus, a

simple example problem is presented to demonstrate these effects. A square, two-dimensional (2-D), 

plane strain composite section was loaded in uniaxial tension, and failure was allowed to evolve within 

microscale RUCs that were linked to the macroscale with FEAMAC. The overall SERR of the model was 

calculated when the dimensions of the RUC were scaled such that the length of the sides of the RUC are 

equal to the characteristic length of the macroscale finite element, and the SERR of the multiscale model 

was also calculated when the dimensions of the RUC were fixed to a realistic value. The calculations 

were performed for numerous macroscale FEM mesh densities and shapes (including square and 

triangular).

2.0 Modeling Progressive Failure at the Microscale

A micromechanical analysis technique, coined the method of cells (MOC), was developed by Aboudi 

(1991); subsequently, Paley and Aboudi (1992) expanded MOC into the generalized method of cells 

(GMC), and later Aboudi et al. (2001) further increased the robustness of this method with the high 

fidelity generalized method of cells (HFGMC). These methods provide semi-closed form solutions for 

determining global anisotropic composite properties in terms of the constituent materials, as well as the 

full three dimensional (3-D) stresses and strains in each of the constituent subcells. 

2.1 The Generalized Method of Cells

It is assumed that a continuously reinforced composite microstructure can be represented as a 

collection of doubly-periodic RUCs containing a general number of constituents, as shown in Figure 1. 

The RUC is then discretized into Nβ×Nγ rectangular subcells, as exhibited in Figure 2. Each of these 

subcells is occupied by one of the constituent materials of the composite. The number of subcells and the 

number of materials are completely general. 

Figure 1.—Representation of the 

doubly-periodic microstructure 

of a composite material.
NASA/TM—2012-217650 3

Figure 2.—Discretization of a doubly-periodic RUC.

Implementation of a Smeared Crack Band Model.

Displacement and traction continuity are enforced in an average, integral sense at the subcell 

interfaces of a discretized RUC. These continuity conditions are used to formulate a set of semi-analytical 

linear algebraic equations that are solved for the local subcell strains in terms of globally applied strains 

or stresses. Subsequently, local constitutive laws can be utilized to obtain the local subcell stresses, and a 

constitutive law for the effective, homogenized composite can be formulated. GMC and HFGMC were 

later reformulated reducing the total number of unknowns required to define the problem by Pindera and 

Bednarcyk (1997) and Aboudi et al. (2012), respectively.

2.2 Implementation of a Smeared Crack Band Model

When materials enter the post-peak strain softening regime (i.e., failure), stress is related to strain 

through a negative tangent stiffness. If the tangent stiffness tensor looses positive definiteness, failure 

localizes to the smallest length scale in the continuum problem (Bažant and Cedolin, 1991). In FEM this 

is a single element, and in GMC this is a single subcell. Thus, the post-peak strain softening energy is 

dissipated over the volume of the discretization entity it localizes to. Since a stress-strain relationship 

prescribes the energy density dissipated during the failure process, the total amount of energy dissipated 

in the local volume decreases if the size of that local volume is reduced, and in the limit the structure 

would theoretically require zero energy to fail. This pathological mesh dependence of numerical 

simulations utilizing constitutive laws containing a negative tangent stiffness that relates stresses to 

strains is well documented in the literature (Bažant and Cedolin, 1979; Pietruszczak and Mroz, 1981; 

deBorst, 1987).

One way to overcome this mesh dependence is the use of non-local constitutive laws, which utilize 

the strain gradients to introduce the influence of surrounding elements into the stress-strain behavior of a 

failing element, preventing localization (Eringen, 1996; Bažant, 1994). A more tractable method involves 

judiciously scaling the post-peak softening slope of the stress-strain curve. Then, the energy density 

dissipated via failure is a function of the characteristic length within the discretized continuum, and the 

strain energy release rate (SERR) (energy required to create a unit area of new crack surface) is preserved 

despite changes in the size of the mesh. This type of constitutive law is referred to as a crack band, or 

smeared crack, model (Bažant and Oh, 1983; de Borst and Nauta, 1985).

A variation of the crack band model was implemented in HFGMC by Pineda et al. (2012), and was 

also implemented in GMC for this work. Crack band initiation is determined using a maximum principal 

stress criterion.
NASA/TM—2012-217650 4

1, 1 ( )

( ) 1

( )

1 = ≥ βγ

βγ

βγ

σ

σ

σ

C

(1)

where ( )

1

βγ σ is the maximum principal stress, and (βγ ) σ C is the mode I cohesive strength of the material.

Upon initiation, the crack band is oriented such that the normal to the crack band ( )

1

βγ n is aligned with the 

direction of maximum principle stress. Figure 3 displays a crack band within a discretized continuum, and

a magnified view of a single subcell shows the crack band of width wC oriented perpendicular to ( )

1

βγ n .

The characteristic length of the subcell (βγ )

C

l is taken as the length of a line running across the subcell, 

oriented in the direction of ( )

1

βγ n .

The behavior of the crack band is controlled by a traction t versus separation [u] law, shown in 

Figure 4(a). The traction-separation law prescribes that the SERR of the subcell (given by the area under 

the traction separation curve) is equal to the fracture toughness, or critical SERR (GC), of the material 

when the traction reaches zero. Since no actual discontinuity is present within the GMC subcell, the 

subcell behaves as a continuum and obeys a stress versus strain constitutive law, displayed in Figure 4(b).

The dissipated strain energy density WF is the area under the stress-strain curve. Additionally, the stress 

state in a subcell is uniform; therefore, the subcell stresses and tractions on the crack band are equivalent.

To ensure mesh objectivity, the negative, post-peak softening slope of the stress-strain curve ET is 

adjusted so that, for subcell βγ which contains the crack band,

( ) (βγ ) (βγ ) (βγ ) (βγ )

C F C GIC l W l = (2)

Figure 3.—(Left) Crack band embedded in a discretized continuum. 

(Right) Magnified subcell displays crack band orientation within 

subcell, as well as characteristic length of subcell.

(a) (b)

Figure 4.—(a) Traction separation law governing behavior of crack 

band. (b). Stress-strain behavior of material.
NASA/TM—2012-217650 5

where (βγ )

C

l is the characteristic length of the subcell, defined in Figure 3, (βγ ) WF is the total strain energy 

density, which is dependent on (βγ )

C

l ,

that is dissipated in the subcell when the stress reaches zero, and 

(βγ ) GIC is the fracture toughness of the material occupying the subcell (material property). Equation (2)

ensures objectivity with respect to changes in the subcell mesh by maintaining that the total SERR, upon 

complete failure (i.e., completely diminished stress), in the subcell is fixed and independent of the size of 

the subcell. For complete details on the implementation of the crack band model, please refer to Pineda et 

al. (2012).

3.0 FEAMAC Multiscale Framework

A multiscale framework is utilized to simulate the deformational response of fiber-reinforced 

composite structures by modeling the fiber-matrix architecture as an RUC at the microscale using GMC

and coupling the microscale to the lamina/laminate level (macrsocale) FEM model. A synergistic 

approach is employed which executes concurrent multiscaling in time, but sequential in length (Sullivan 

and Arnold, 2011). The commercial finite element software, Abaqus 6.10-1 (Simulia Corp., 2008), is used 

as the FEM platform, and the MAC/GMC core micromechanics software (Bednarcyk and Arnold, 2002a; 

Bednarcyk and Arnold, 2002b) is used to perform microscale calculations. The scales are linked using the 

FEAMAC software implementation (Bednarcyk and Arnold, 2006), which utilizes various 

Abaqus/Standard user subroutines. Subsequently, a wrapper was developed (FEAMAC/Explicit) which 

achieves the same multiscaling tasks as FEAMAC, except Abaqus/Explicit is used as the FEM platform 

(Pineda et al., 2009). MAC/GMC, FEAMAC, and FEAMAC/Explicit are available in the openly 

distributed ImMAC suite of micromechanics codes, developed by GRC. The number of subcells and 

constituents within an RUC at the microscale is completely general, and as many, or as few, as needed 

may be used to accurately characterize the micro-architecture of the composite. Moreover, any 

constitutive model can be used in the subcells of the RUC to rerpresent the mechanical or thermal 

response of the constituents.

A schematic displaying a typical multiscale model using FEAMAC is displayed in Figure 5. The 

integration point strains are applied to the RUC and the local subcell fields are determined using GMC 

(this process is referred to as localizaion). If the subcell material behavior is nonlinear, the local stresses 

and strains are used to calculate the local stiffnesses, plastic strains, thermal strains, and/or state variables

via the local constitutive law. The RUC is then homogenized and the global stiffnesses, plastic strains, 

Figure 5.—Diagram showing coupling of macroscale FEM and microscale GMC models.
NASA/TM—2012-217650 6

thermal strains, and/or state variables are computed. The global stresses at the integration point are then 

calculated using these global, homogenized fields, and the effects of any nonlinear subscale phenomena 

are introduced into the macroscale through changes in the integration point stress state. The global 

stresses, material jacobian, and updated state variables are then supplied to the Abaqus user material 

UMAT subroutine. For complete details on the FEAMAC implementation, the reader is referred to 

Bednarcyk and Arnold (2006).

4.0 Effect of Subscale Size on Global Energy Dissipation

A mesh objective continuum damage mechanics (CDM) model for post-peak strain softening

implemented within the GMC micromechanics framework was presented in Section 2.2. Using such a 

model ensures that the solutions obtained are insensitive to changes in the subcell mesh (Pineda et al., 

2012). Nonetheless, this does not guarantee mesh objectivity if this model is linked to a macroscale model 

via a multiscale framework. To ensure mesh objectivity the SERR of all elements must be preserved upon 

changes in the macroscale mesh.

Typically, an RUC is considered dimensionless as the periodic boundary conditions are meant to 

simulate an infinite array of the RUC. However, the domain of the integration point that the microscale 

model is coupled to has some finite dimensions. As energy is dissipated at the microscale within the 

RUC, the macroscale also exhibits energy dissipation governed by the stress-strain response of the 

homogenized RUC, and this energy is dissipated over the volume of the integration point domain.

The dissipated energy density at the microscale RUC WD can be calculated as the total energy density 

minus the elastic strain energy density.

RUC RUC

ij

RUC

ij

RUC

ij

RUC

RUC ij

RUC

D d dV

V

W ∫ ∫ 



 



 = σ ε − σ ε

2

1 1

(3)

where RUC

ij ε and RUC σ ij are the global, homogenized RUC strain and stress tensor, respectively, and RUC V

is the volume of the RUC. Thus, the SERR RUC G can be calculated as the total dissipated energy density 

multiplied by the characteristic length of the RUC.

RUC

D

RUC

C

RUC G = l W (4)

If the volume of the RUC remains constant for all simulations, then the SERR also remains constant.

Noting that, in a multiscale model, the integration point stress and strains are equivalent to the global 

stresses and strains of the homogenized RUC, then the dissipated strain energy density is automatically 

preserved across the scales.

RUC

D

I P WD =W . . (5)

where I.P. WD is the dissipated energy density of the integration point. The SERR of the integration point 

domain I.P. G can be calculated by multiplying Equation (3) by the characteristic length of the integration 

point domain I.P.

Cl .

RUC

D

I P

C

I P G l W . . . . = (6)

In multiscale simulations that keep the volume of the RUC fixed, the dissipated strain energy density of 

the RUC RUC WD is constant, and thus, the SERR is directly proportional to the characteristic length of the 

integration point domain to which the failure localizes. Consequently as the mesh is refined, the SERR in 
NASA/TM—2012-217650 7

the localized domains, and thus the entire structure, changes. The results will never converge and are 

pathologically mesh dependent.

To enact mesh objectivity in the multiscale model, the RUC must preserve the characteristic length of 

the integration point domain.

I.P.

C

RUC

Cl = l (7)

If the characteristic length is preserved, then Equation (6) becomes

I P RUC G = G . . (8)

and the SERR of the RUC is preserved within the integration point domain after homogenization. Since 

the local SERR release rate in the failing subcells is governed by the critical fracture toughness (a fixed 

material property), the global SERR of the RUC is independent of size. By extension, the SERR of the 

integration point domain is also insensitive to changes in size. 

The challenge in achieving mesh objective multiscale analyses lies in correctly defining the global 

characteristic length RUC

Cl , or I.P.Cl . Locally within the RUC, the characteristic length can be defined for a 

subcell by the orientation of the crack band in the subcell (see Section 2.2). Conversely, in the

homogenized RUC, there is no single crack band angle to define the characteristic length. This can be 

assuaged by defining the RUC such that the dimensions are exactly the same as the dimensions of the 

integration point domain. Thus, every length that can be defined in the integration point domain is defined 

equivalently in the RUC, and the characteristic length is automatically preserved across the scales, even if 

it is unknown.

Unfortunately, in the doubly-periodic GMC micromechanics framework, the shapes of the RUCs are 

restricted to rectangles; however, 2-D finite elements may take the shape of any general quadrilateral. The 

following sections investigate the effects of scaling the size of the RUC to correspond to the volume of 

the FEM integration point domain on the total SERR of the multiscale model. To demonstrate mesh 

objectivity, the total SERR of the entire model must be preserved as the FEM mesh is changed.

4.1 Numerical Model

A simple example problem is used to demonstrate the effects of scaling the dimensions of the 

microscale RUC appropriately on the overall SERR of the multiscale model. Additionally, the effects of 

using incompatible finite elements and RUC shapes are investigated with the same example problem.

Figure 6 exhibits the multiscale example problem. At the macroscale a 1- by 1-mm square section is 

modeled using a uniform mesh consisting of square, 2-D, plane strain, reduced integration, CPE4R 

elements with dimensions le x le. A global x-y-z coordinate system is employed at the macroscale. The 

microscale domain is modeled using a doubly-periodic, 7 subcells x 7 subcells GMC RUC with overall 

dimensions L x L. To maintain the proper aspect ratio, for all RUC sizes, to represent a circular fiber the 

height and length of the RUC must be equal. The RUC consists of 13 fiber subcells (colored blue) and 36 

matrix subcells (colored green), and is meant to represent a square packing architecture. A local x1-x2-x3

coordinate system is used where x1 is the axial, fiber direction, and x2/x3 are the transverse directions. The 

x1/z-directions are coincident, as are the x2/x-directions, and the x3/y-directions.

The elastic properties of the fiber and matrix constituent subcells were chosen so that the overall 

elastic response of the homogenized RUC represented the elastic behavior of an IM7/8552 lamina. The 

transversely elastic properties of IM7/8552 are presented in Table 1 (Camanho et al., 2007), where the x-y 

plane in Figure 6 is the plane of isotropy and Ezz is the axial (fiber direction) Young’s modulus, Exx is the

transverse stiffness, νzx is the axial Poisson’s ratio, and Gzx is the axial shear modulus. The matrix subcells 

in the GMC RUC initially behave as an elastic, isotropic material, and the fibers subcells are initially 

elastic and transversely isotropic. The initial, elastic properties of the microscale constituents are given in 
NASA/TM—2012-217650 8

Table 2, where Em is the matrix Young’s modulus, ν

m is the matrix Poisson’s ratio, f E11 is the axial fiber

Young’s modulus in the x1-direction (see Figure 6), f E22 is the transverse stiffness of the fiber, f ν12 and 

f

23 ν are the axial and transverse Poisson’s ratios, and f G12 is the axial shear stiffness. A fiber volume 

fraction Vf

of 59.1 percent was used (Camanho et al., 2003). The transverse and axial shear moduli of the 

fiber were obtained from Goldberg and Gilat (2003). The elastic properties of the matrix and remaining 

unknown elastic properties of the fiber were calibrated such that the elastic stiffness of the homogenized 

RUC correlated to the elastic properties of the IM7/8552 lamina.

Figure 6.—Multiscale simulation utilizing square elements at the macroscale

TABLE 1.—ELASTIC PROPERTIES OF 

IM7/8552 LAMINA (CAMANHO ET AL., 2007)

Ezz, GPa................................................................. 171.4 

Exx, GPa................................................................... 9.08 

νzx ............................................................................ 0.32

Gzx, GPa .................................................................. 5.29

TABLE 2.—ELASTIC PROPERTIES OF IM7 CARBON FIBER AND 

8552 EPOXY MATRIX CONSTITUENTS USED 

IN MICROSCALE GMC RUC

8552 Matrix Properties

Em, GPa ................................................................................. 4.97

ΝM.......................................................................................... 0.36

IM7 Fiber Properties

f E11 , GPa ............................................................................. 286.5 

f E22 , GPa............................................................................... 12.4 

f ν 12 ........................................................................................ 0.29

f ν 23 ........................................................................................ 0.29

f G12 , GPa ............................................................................... 20.0
NASA/TM—2012-217650 9

The bottom edge of the FEM domain was pinned in the y-direction, and the left edge was pinned in 

the x-direction. A uniform displacement was applied to the right edge of the FEM model to simulate 

transverse loading of a 1- by 1-mm composite section that stretches infinitely in the global, axial zdirection (plane strain conditions). Due to the high degree of nonlinearity associated with failure, the 

implicit, dynamic FEM solver available in Abaqus/Standard was utilized to solve the macroscale FEM 

problem using the keyword *DYNAMIC, with the option APPLICATION = QUASI-STATIC (Simulia 

Corp., 2008). A nominal density of 1.57 kg/m3 was used in all the finite elements. The total displacement 

(ranging from 0.0175 to 0.05 mm), was applied over a time of 1000 sec.

Failure is introduced into the matrix subcells of the microscale RUC, and the effects are transferred to 

the macroscale through an overall reduction of the stiffness of the RUC (due to a local reduction in the 

stiffness of the failing matrix subcells). The crack band model presented in Section 2.2 is used to govern 

failure at the microscale. Only the yellow and red elements in Figure 6 were modeled with multiscale 

FEAMAC elements. The green elements were modeled using the Abaqus elastic material model and the 

properties listed in Table 1. Free edge effects and the accuracy of the local fields surrounding the failure 

path front are influenced by the mesh refinement, to some degree, which may affect the failure path in the 

simulation. Therefore, failure was restricted to the yellow and red elements to simplify the problem as 

much as possible and try to ensure that the energy dissipated in the model is a result of the multiscale 

methodology, not other influences. Moreover, the matrix subcells of the RUC located in the red element 

were given a slightly lower cohesive strength than the matrix of the RUCs within the yellow elements.

This was imposed to imitate an imperfection and facilitate the localization of initial failure into a single 

element. Without this imperfection, all of the FEAMAC elements would fail simultaneously because the 

stress and strain fields are uniform throughout the FEM mesh. After failure initiates in the RUC located in 

the red element, failure progresses into the RUCs located within the adjacent yellow elements until the 

failure reaches the free edge of the FEM mesh.

The local mode I cohesive strength (βγ ) σ C , and the mode I fracture toughness (βγ ) GIC of the matrix 

subcells were calibrated using a macroscale mesh of 31 elements x 31 elements. The dimensions of the 

microscale RUC were chosen to equal the dimension of the finite element; i.e., L = l

e = 0.0323 mm. The 

local parameters were calibrated so that the global ultimate stress and total SERR were close to 60.3 MPa

and 0.2774 kJ/m2

, respectively, which are the transverse strength and mode I, transverse fracture 

toughness for IM7/8552 reported by Camanho et al. (2007). The mode I cohesive strength, and mode I 

fracture toughness of the matrix that were used in the simulations are displayed in Table 3. It should be 

noted, that the critical strain of the matrix subcells of the RUCs located in the yellow elements were given 

a value that was 5 percent higher than that used in the red element.

TABLE 3.—CRACK BAND FAILURE PARAMETERS 

USED IN 8552 MATRIX SUBCELLS

Mode I cohesive strength (βγ ) σ C (Weak) .............56.7 MPa (54.0 MPa)

Mode I fracture toughness (βγ ) GIC ........................................ 1.594 kJ/m2

Simulations were executed using various square meshes, and the total, macroscopic SERR was 

calculated upon complete diminishment of the load carrying capability. Two scenarios were investigated.

The first involved scaling the dimensions of the microscale RUC, so that

I.P.

C L = l (9)

while retaining Vf= 59.1 percent. Since reduced integration elements are used, the characteristic length of 

an element e

Cl , and of the integration point domain in an element I.P.

Cl , are equivalent.
NASA/TM—2012-217650 10

e

C

I P

Cl = l . . (10)

The characteristic length of a 2-D element is calculated as the square root of the area of the element Ae

(Simulia Corp., 2008).

e e lC = A (11)

Thus, for a square, reduced integration element, the characteristic length of the integration point domain 

within an element is equal to the length of one of the sides of the element.

I P e

Cl = l . . (12)

In the second scenario, the dimension of the RUC was kept fixed at L = 5.98 μm (this was calculated 

assuming a fiber diameter of 5 μm, and retaining a 59.1 percent fiber volume fraction), while the FEM 

mesh was refined.

To determine an upper bound on the error introduced by having incompatible element and RUC 

shapes, the effects of using 2-D, plane strain, triangular CPE3 elements on the total SERR of the model 

was also investigated. Figure 7 displays a typical macroscale mesh utilizing triangular elements. The 

dimensions of the RUC were still scaled according to Equation (9); however, the definition of the 

characteristic length of the integration point domain is no longer given by Equation (12). Utilizing 

Equation (11), I.P.

Cl becomes

I P e

Cl l

2

. . 2 = (13)

for triangular elements. Additionally, the same simulations were conducted with fixed RUC dimension 

L = 5.98 μm.

Figure 7.—Multiscale simulation utilizing triangular elements at the macroscale.
NASA/TM—2012-217650 11

4.2 Results

The total SERR of the model G was calculated using

rightedge

rightedge

rightedge

rightedge

rightedge

rightedge

U

L

RF

dU

L

RF

G ∑ ∫ = ∑ − 2

1 (14)

where ∑ rightedge RF is the sum of the reaction forces at all the nodes on the right edge of the model, 

Lrightedge is the length of the right edge (1 mm for all simulations), and Urightedge is the uniform displacement 

applied to the right edge of the model. Figure 8 displays G as a function of the characteristic length I.P.

Cl

for the various meshes. The results for the square meshes, shown in Figure 6, are indicated with square 

markers; whereas, the results for the triangular meshes (see Figure 7) are indicated with triangular 

markers. Filled markers indicate that the dimensions of microscale RUCs were scaled to satisfy 

Equation (9); whereas empty markers indicated that the RUC dimensions were fixed, L = 5.98 μm.

The element length that is closest to the realistic size of the RUC (L = 5.98 μm) is e

Cl = 5.99 μm 

(167 elements x 167 elements), and is only 0.17 percent larger than the actual RUC length. Thus, the 

scaled solution obtained with I.P.

Cl = 5.99 μm is considered to be the most accurate solution. For this mesh 

G was calculated to be 0.227 kJ/m2

; this value will be used as the baseline SERR for comparison with the 

other simulations. The most erroneous G was calculated, from the solution using the coarsest mesh

I.P.

Cl = 32.3 μm, to be 0.263 kJ/m2

. Therefore, for the worst solution, G was 16 percent higher than the 

most accurate solution; although, the element length was 438.7 percent larger. As the mesh was refined, G

approached 0.227 kJ/m2. The maximum G, minimum G and percentage difference between the two Gs 

and corresponding I.P.

Cl s are presented in Table 4 for all the square and triangular mesh simulations.

Figure 8.—Total SERR of multiscale models as a function of 

characteristic length.
NASA/TM—2012-217650 12

TABLE 4.—MAXIMUM AND MINIMUM G, AS WELL AS CORRESPONDING I.P.

Cl , 

AND PERCENTAGE DIFFERENCE BETWEEN GS AND I.P.

Cl

Quad. - scaled Quad. - unscaled Tri. - scaled Tri. - unscaled

I.P.

Cl ,

μm

G,

kJ/m2

I.P.

Cl ,

μm

G

kJ/m2

I.P.

Cl ,

μm

G,

kJ/m2

I.P.

Cl ,

μm

G,

kJ/m2

Max G 32.3 0.263 32.3 0.702 33.7 0.397 33.7 0.866

Min G 5.99 0.227 5.99 0.230 5.94 0.276 6.15 0.281

% Diff. 438.7 16.0 438.7 204.8 467.3 44.0 448.0 208.2

If the dimensions of the RUC were not scaled, the error observed in the SERR was very significant.

The coarsest mesh exhibited G = 0.702 kJ/m2 which is an error of 209.6 percent from the baseline.

Although, when the characteristic length of the integration point domain and the size of the RUC are 

comparable (i.e., I.P.

Cl = 5.99 μm), G = 0.230 kJ/m2

, which is only difference of 1.6 percent from the 

baseline. However, the SERR calculated for the unscaled simulations dropped suddenly when the element 

size approached the size of the RUC. For example, the average G for the next three largest meshes 

( I.P.

Cl = 6.14, 7.41, and 8.70 μm) is 0.328 kJ/m2

. The SERR for the coarsest mesh is 204.8 percent larger 

than the finest mesh if the RUC dimensions are not scaled.

Clearly, scaling the size of the RUC such that it corresponds to the dimensions of the finite element 

minimizes the error. The variation that is observed in simulations with scaling, as the mesh is refined, can 

be attributed to local effects in the macroscopic FEM model such as an increase in the stress 

concentration ahead of the failure path front, as the mesh is refined.

The simulations incorporating triangular meshes exhibited an overall larger degree of error (when 

compared to the baseline) in the SERR. This is to be expected, since the shape the element and the shape 

of the RUC are incompatible. The characteristic length, calculated using Equation (13), that was closest to 

the assumed actual RUC length, is 5.94 μm. When the RUC dimensions were scaled, G was calculated to 

be 0.276 kJ/m2 for this characteristic element length, an error of 21.6 percent from the baseline calculated 

using square elements. The largest G calculated (when RUC scaling was enforced), corresponding to the 

coarsest mesh with I.P.

Cl = 33.7 μm, was 0.397 kJ/m2

. The difference between the minimum G and 

maximum G is 44.0 percent, nearly three times the difference seen in simulations with square elements.

For the unscaled simulations utilizing a triangular mesh, the largest G calculated was 0.866 kJ/m2

,

corresponding to the coarsest mesh ( I.P.

Cl = 33.7 μm), and the smallest G calculated was 0.281 kJ/m2

, 

corresponding to I.P.

Cl = 6.15 μm. Note that the finest mesh ( I.P.Cl = 5.94 μm) produced an SERR of 

0.294 kJ/m2 which is slightly higher than the SERR obtained from the slightly coarser mesh. The 

minimum error was 23.9 percent; however the difference between the maximum G and the minimum G

was 208.2 percent. The mismatch between the element shape and the RUC shape introduces some error 

when compared to the baseline square solution; thus, the use of triangular elements in an FEAMAC 

simulation should be avoided. However, the error introduced through refining the mesh can be minimized 

by scaling the RUC dimension by the characteristic length of the element.

5.0 Conclusions

A smeared crack band model was implemented within the GMC micromechanics framework. The 

crack band model utilizes the characteristic length of the subcells in the micromechanics RUC to preserve 

the SERR of each subcell, ensuring that it is equal to the critical fracture toughness of the material and is 

insensitive to changes in the subcell mesh. This yields mesh objective microscale failure.

A platform for linking macroscale finite element analysis to microscale GMC computations, the 

FEAMAC framework developed at GRC, is presented. This framework utilizes standard Abaqus and 

MAC/GMC input decks and data cards to setup a multiscale problem. Communication between the 
NASA/TM—2012-217650 13

MAC/GMC core engine, also developed at GRC, and Abaqus FEM commercial software is facilitated 

through an Abaqus user material UMAT subroutine.

Use of the crack band model at the microscale ensures the microscale behavior is objective with 

respect to changes in the GMC subcell mesh; however, it does not guarantee that the macroscopic 

solution will be objective with respect to the finite element mesh in a multiscale simulation. A scaling 

technique was developed for preserving the inherent length scale of the macroscale as the analysis 

proceeds into the microscale. The scaling procedure utilizes the characteristic length of the integration 

point domain within an element to scale the dimensions of the RUC such that the volume of the 

integration point domain and the RUC are equivalent. Plane strain, 2-D, uniaxial tension simulations were 

performed using meshes with square, triangular, and rectangular elements with, and without, scaling the 

RUCs for numerous degrees of mesh refinement. The overall SERR, compared to a baseline model, of 

each simulation is used to gauge the accuracy of the models. The technique works well when the elements 

are nearly square in shape, and alleviates a large degree of error observed when the scaling is not 

performed. Solutions using triangular elements exhibited significant error; however, the results were 

fairly insensitive to mesh refinement. However, changing the aspect ratio of rectangular elements 

significantly affected the accuracy of the solutions, and scaling the RUC by the characteristic length of 

the elements did not appropriately adjust the microscale problem to preserve the SERR as the mesh was 

changed.

References

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

Approach (in press), Elsevier, 2012.

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

Amsterdam, Netherlands, 1991.

Aboudi, J., Pindera, M.–J., and Arnold, S. M., “Linear Thermoelastic Higher-Order Theory for Periodic 

Multiphase Materials,” J. Appl. Mech., vol. 68, pp. 697-707, 2001.

Bažant, Z. P., “Can Multiscale-Multiphysics Methods Predict Softening Damage and Structural Failure?” 

Mech. Amer. Acad. Mech., vol. 36 (5-6), pp. 5-12, 2007.

Bažant, Z. P., and Oh, B. H., “Crack Band Theory for Fracture of Concrete,” Mater. Struct., vol. 16, pp. 

155-177, 1983.

Bažant, Z.P., and Cedolin, L., “Blunt Crack Band Propagation in Finite Element Analysis,” J. Eng. Mech. 

Div.-ASCE, vol. 105, pp. 297-315, 1979.

Bažant, Z.P., and Cedolin, L., Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories, 

Oxford University Press, New York, NY, 1991.

Bednarcyk, B. A., and Arnold, S. M., “A Framework for Performing Multiscale Stochastic Progressive 

Failure Analysis of Composite Structures,” Proceedings of the 2006 Abaqus User’s Conference, May 

23-25, Boston, MA, 2006.

Bednarcyk, B. A., and Arnold, S. M., “MAC/GMC 4.0 User’s Manual – Keywords Manual,” 

NASA/TM—2002-212077/VOL2, 2002a.

Bednarcyk, B. A., and Arnold, S. M., “MAC/GMC 4.0 User’s Manual – Example Problems Manual,” 

NASA/TM—2002-212077/VOL3, 2002b.

Camanho, P. P., Maimí, P., and Dávila, C. G., “Prediction of Size Effects in Notched Laminates Using 

Continuum Damage Mechanics,” Compos. Sci. Technol., vol. 67, pp. 2715-2727, 2007.

de Borst, R., and Nauta, P., “Non-orthogonal Cracks in a Smeared Finite Element Model,” Eng. Comput., 

vol. 2, pp. 35-46, 1985.

deBorst, R., “Computation of Post-bifurcation and Post-failure Behavior of Strain-softening Solids,” 

Comput. Struct., vol. 25, pp. 211-224, 1987.

Goldberg, R. K., and Gilat, A., “Experimental and Computational Characterization of the High Strain 

Rate Tensile Response of Polymer Matrix Composites,” Composite Materials: Testing and Design 
NASA/TM—2012-217650 14

Fourteenth Volume, ASTM STP 1436, C. E. Bakis, Ed., ASTM International, West Conshohocken, 

PA, 2003.

Paley, M., and Aboudi, J., “Micromechanical Analysis of Composites by the Generalized Cells Model,” 

Mech. Mater., vol. 14, pp. 127-139, 1992.

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

J. Numer. Methods Eng., vol. 17, pp. 327-334, 1981.

Pindera, 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.

Pineda, E. J., Waas, A. M., Bednarcyk, B. A. and Arnold, S. M., “Implementation of a Smeared Crack 

Band Model in a Micromechanics Framework,” NASA/TM (in-press), 2012.

Pineda, E. J., Waas, A. M., Bednarcyk, B. A., Arnold, S. M., and Collier, C. S., “Multiscale Failure 

Analysis of Laminated Composite Panels Subjected to Blast Loading Using FEAMAC/Explicit,” 

NASA/TM—2009-215813, 2009.

Simulia Corp., Abaqus User’s Manual, Vol. 1-3, Version 6.10-1, Providence, RI, 2008.

Sullivan, R. W. and Arnold, S. M., “An Annotative Review of Multiscale Modeling and its Application to 

Scales Inherent in the Field of ICME,” Models, Databases, and Simulation Tools Needed for the 

Realization of Integrated Computational Material Engineering, S. M. Arnold and T. T. Wong, Eds., 

ASTM International, Materials Park, OH, pp. 6-23, 2011.

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

Effects of Subscale Size and Shape on Global Energy Dissipation in a Multiscale Model of a 

Fiber-Reinforced Composite Exhibiting Post-Peak Strain Softening Using Abaqus and 

FEAMAC 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Pineda, Evan, J.; Bednarcyk, Brett, A.; 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-18328 

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

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 24, 39, and 61 

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

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

13. SUPPLEMENTARY NOTES

14. ABSTRACT

A mesh objective crack band model is implemented in the generalized method of cells (GMC) micromechanics model to predict failure of a 

composite repeating unit cell (RUC). The micromechanics calculations are achieved using the MAC/GMC core engine within the ImMAC 

suite of micromechanics codes, developed at the NASA Glenn Research Center. The microscale RUC is linked to a macroscale 

Abaqus/Standard finite element model using the FEAMAC multiscale framework (included in the ImMAC suite). The effects of the 

relationship between the characteristic length of the finite element and the size of the microscale RUC on the total energy dissipation of the 

multiscale model are investigated. A simple 2-D composite square subjected to uniaxial tension is used to demonstrate the effects of scaling 

the dimensions of the RUC such that the length of the sides of the RUC are equal to the characteristic length of the finite element. These 

results are compared to simulations where the size of the RUC is fixed, independent of the element size. Simulations are carried out for a 

variety of mesh densities and element shapes, including square and triangular. Results indicate that a consistent size and shape must be used 

to yield preserve energy dissipation across the scales. 

15. SUBJECT TERMS

Fiber composites; Polymers matrix composites; Carbon fiber reinforced plastics; Computational mechanics; Damage; Fracture 

mechanics 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

22 

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