Vipul Ranatunga

Miami University, Middletown, Ohio

Brett A. Bednarcyk and Steven M. Arnold

Glenn Research Center, Cleveland, Ohio

Modeling Progressive Damage Using Local

Displacement Discontinuities Within the

FEAMAC Multiscale Modeling Framework

NASA/TM—2010-216825

October 2010

AIAA–2010–2619
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
Vipul Ranatunga

Miami University, Middletown, Ohio

Brett A. Bednarcyk and Steven M. Arnold

Glenn Research Center, Cleveland, Ohio

Modeling Progressive Damage Using Local

Displacement Discontinuities Within the

FEAMAC Multiscale Modeling Framework

NASA/TM—2010-216825

October 2010

AIAA–2010–2619

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

51st Structures, Structural Dynamics, and Materials Conference 

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

Orlando, Florida, April 12–15, 2010
Acknowledgments

The fi rst author would like to recognize the support received through the NASA Summer Faculty Fellowship Program, 

supported by the Exploration Technology Development Program (ETDP) Advanced Composites Technologies (ACT)

Project and the software and hardware resources utilized for this research from the

Research Computing Support Group at Miami University.

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://gltrs.grc.nasa.gov

Level of Review: This material has been technically reviewed by technical management. 
NASA/TM—2010-216825 1 

Modeling Progressive Damage Using Local Displacement 

Discontinuities Within the FEAMAC Multiscale 

Modeling Framework 

Vipul Ranatunga 

Miami University 

Middletown, Ohio 45042 

Brett A. Bednarcyk and Steven M. Arnold 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Abstract 

A method for performing progressive damage modeling in composite materials and structures based on 

continuum level interfacial displacement discontinuities is presented. The proposed method enables the exponential 

evolution of the interfacial compliance, resulting in unloading of the tractions at the interface after delamination or 

failure occurs. In this paper, the proposed continuum displacement discontinuity model has been used to simulate 

failure within both isotropic and orthotropic materials efficiently and to explore the possibility of predicting the 

crack path, therein. Simulation results obtained from Mode-I and Mode-II fracture compare the proposed approach 

with the cohesive element approach and Virtual Crack Closure Techniques (VCCT) available within the ABAQUS 

(ABAQUS, Inc.) finite element software. Furthermore, an eccentrically loaded 3-point bend test has been simulated 

with the displacement discontinuity model, and the resulting crack path prediction has been compared with a 

prediction based on the extended finite element model (XFEM) approach. 

I. Introduction 

Carbon fiber-reinforced polymers (CFRPs) are used widely for many applications including aircraft, 

automotive, and civil structures. Some of the appealing properties of CFRPs include high specific strength and 

stiffness, corrosion resistance, fatigue resistance, and suitability for stealth applications. Increased usage of CFRPs 

in the future depends on the development of damage-tolerant composite structures while lowering the cost of 

manufacturing. One of the issues with CFRPs for structural applications is the joining of sections made with fiberreinforced composites. Among other joining methods in current practice, two of the most popular methods are 

adhesive bonding and mechanical fastening where the mechanical fastening is preferred for thicker sections while 

adhesive bonding is used for thinner sections. Delamination is a major type of failure associated with adhesivelybonded composite joints. 

In the past, failure of composites has been investigated extensively from the micromechanical and 

macromechanical points of view. On the micromechanical scale, failure mechanisms are related to the properties of 

the constituent phases, i.e., matrix, reinforcement, and the interface. This paper will investigate the advantages and 

disadvantages of modeling damage propagation in materials utilizing a displacement discontinuity approach. The 

displacement discontinuity approach known as ECI (evolving compliant interface) (Refs. 1 and 2) enables 

exponential evolution of the interfacial compliance, resulting in unloading of the tractions at the interface after 

failure/delaminations occurs. The current study utilizes the implementation of this approach available within 

NASA’s FEAMAC (Ref. 3) multiscale micromechanics analysis framework. Given the fact that the commercially 

available techniques such as VCCT (Ref. 4) and cohesive element (Ref. 5) approaches require a priori knowledge of 

the crack path, and as discussed elsewhere (Ref. 6), the current implementation of these approaches present 

difficulties with numerical convergence, forcing one to deal with a number of parameters to overcome these 

convergence issues. Therefore, in this paper, a computationally-efficient multiscale modeling framework is 

presented, and the propagation of damage is investigated utilizing this framework in the context of finite element 

method. The future intent is to apply this approach to the failure/delaminations of composite structures. 

To begin the investigation, VCCT and cohesive element approaches are used to model self-similar cracks in 

double cantilever beam (DCB) the end-notched flexure (ENF) tests, and the results are compared with FEAMAC 

predictions. Once the material data and the simulation-specific parameters have been identified through these 

preliminary simulations, the focus shifts to an eccentrically-loaded three-point bend test. In the absence of 
NASA/TM—2010-216825 2 

experimental data to validate the FEAMAC prediction for the crack path, the extended finite element method 

(XFEM) (Ref. 7) available in the Abaqus finite element package is utilized for comparison purposes. The study 

concludes with highlighting the advantages and disadvantages of the FEAMAC and XFEM approaches. 

II. Displacement Discontinuity Approach and the Evolving Compliance Interface Model 

The displacement discontinuity approach employed herein follows the basic form proposed by Jones and 

Whittier (Ref. 8), in which a debonded interface is modeled by imposing a discontinuity (i.e., “jump”) in a given 

displacement component across the interface. This basic form can be written, 

I I

j j j   u R     (1) 

Where [uj]

I

 is the discontinuity in displacement component j at interface I, j is the corresponding stress component 

at interface I, and Rj is a proportionality constant that may be thought of as a flexibility for the interface (as it relates 

the interfacial displacement jump to the interfacial stress). This flexible interface model has been applied extensively 

to model fiber-matrix interfacial debonding in composite materials. It was incorporated within the method of cells 

micromechanics model by Aboudi (Ref. 9) and within the generalized method of cell by Sankurathri et al. (Ref. 10). 

Achenbach and Zhu (Ref. 11) added the condition to this flexible interface model that, when in compression, the 

interface must behave as if well-bonded (i.e., Rj = 0). Wilt and Arnold (Ref. 12) added finite bond strength to the 

interface, 

I j

DB  , such that the interface behaves as perfectly bonded until the interfacial stress reaches this 

interfacial strength. Bednarcyk and Arnold (Refs. 1 and 2) further extended the Jones and Whittier (Ref. 8) flexible 

interface model to admit time dependence in the Rj parameter. By allowing this parameter to evolve with time, it is 

possible to enable unloading and redistribution of the interfacial stress upon initiation of debonding or delamination. 

Considering the components normal (n) and tangential (t) to a particular interface separately, the equations for the 

displacement discontinuity components can be written as, 

  ;

I I I I n n n n n DB     u Rt     (2) 

  ;

I I I I t

t t t t DB     u Rt     (3) 

where the time-dependence of the Rj parameter, as well as the finite bond strength of the interface, has been denoted. 

This form of the Jones and Whittier (Ref. 8) flexible interface model has been called the evolving compliant 

interface (ECI) model in order to distinguish it from the previous form that incorporates a constant value for the Rj

parameter. 

Many forms for the time-dependence of the Rj parameter, both implicit and explicit, were examined by 

Bednarcyk and Arnold (Ref. 1). It was determined that explicit exponential time-dependence provides the model 

with the desired ability to unload interfacial stress, permitting redistribution of this stress. Further, the ECI model 

compares well (in a qualitative sense) with several previous interfacial constitutive models that incorporate 

interfacial stress unloading, e.g., Needleman (Ref. 13). The ECI model representation of the debonding parameter 

put forth by Bednarcyk and Arnold (Refs. 1 and 2) is, 

  exp 1 DB

j j DB

j

t t Rt tt

                    

 (4) 

where tDB is the time at which debonding occurs (for the particular interface), and j and Bj are parameters that 

characterize the interfacial time-dependent behavior. The typical effective interfacial behavior represented by the 

ECI model is plotted in Figure 1. Until the interfacial stress reaches the bond strength (DB), the interfacial 

displacement is zero. Then, debonding occurs and the interfacial stress rises and then falls as the interface opens. 
NASA/TM—2010-216825 3 

Also presented in Figure 1 is a comparison of the ECI model with two previous interfacial constitutive models 

that also admit unloading of the interfacial stress. The well-known model of Needleman (Ref. 13) (NI model) is one 

such interfacial representation. While this model does not incorporate finite bond strength, the effective interfacial 

behavior is qualitatively similar to that of the ECI model. The fact that the NI model equations are explicitly 

nonlinear in their relation between interfacial tractions and interfacial displacements, can add difficulty to the global 

solution when incorporated within micromechanics models. Further, the effective interfacial response of the 

statistical interfacial failure (SIF) model, developed by Robertson and Mall (Ref. 14), is plotted in Figure 1 as well. 

Its qualitative similarity to the ECI model and NI model is obvious. Since the SIF model, like the NI model, is 

nonlinear, Robertson and Mall (Ref. 14) incorporated a linear approximation of the SIF model (also plotted in 

Fig. 1) within their employed micromechanics model (the method of cells). 

The displacement discontinuity approach with evolving compliance interface (ECI) model available within 

NASA’s FEAMAC micromechanics analysis framework is primarily developed for modeling debonding between 

fiber and matrix materials. The FEAMAC framework and the ECI model are employed in the present investigation 

to examine crack propagation in the context of finite element simulations. Detailed comparisons of the FEAMAC 

predictions against finite element-based VCCT and cohesive element predictions under Mode-I and Mode-II are 

discussed and the advantages and disadvantages of each technique are discussed in the subsequent sections. 

Prediction of the crack path by FEAMAC is compared with XFEM predictions in the absence of experimental data. 

III. Results 

A. Simulations With a DCB Specimen 

The double cantilever beam test has widely been used for obtaining Mode-I interlaminar fracture toughness in 

composite materials. The test uses a geometry shown in Figure 2 with an initial delamination crack and the 

delamination is forced to grow by pulling apart the two plies of the specimen. To begin the comparison between the 

FEAMAC, VCCT, and the cohesive element approaches, a transversely-isotropic material and the geometry shown 

in Figure 2 is considered. The selected beam geometry has a length (L) of 4 in., width (2W) of 0.12 in. and thickness 

of 0.3 in. as used by Song et al. (Ref. 15). The beam includes two sub-laminates, each with a thickness of 0.06 in. 

and an initial crack length (CL) of 1.15 in. between the two sub-laminates. Material properties of AS4/3501-6 

(Ref. 15) given in Table 1 are used in the plane-strain finite element models, considering unidirectional fibers 

aligned with the length direction (Y) of the specimen. 

Figure 1.—Comparison of the ECI model with the Needleman (Ref. 13) 

interface (NI) model and the statistical interfacial failure (SIF) model 

developed by Robertson and Mall (Ref. 14). 

0

40

80

120

160

200

0.000 0.002 0.004 0.006 0.008 0.010

Interfacial Stress (MPa)

Normalized Interfacial Displacement

NI Model

SIF Model

SIF Model Linear 

Approximation

ECI Model

DB (SIF Approximation)

σDB (SIF, ECI)
NASA/TM—2010-216825 4 

Figure 2.—Geometry of the DCB specimen with applied boundary conditions. 

TABLE 1.—MECHANICAL PROPERTIES OF AS4/3501-6 INCLUDING PARAMETERS USED WITH 

TRACTION-SEPARATION BASED COHESIVE ELEMENT SIMULATIONS (REF. 15) 

E1/ksi E2/ksi E3/ksi G12/ksi G13/ksi G23/ksi 

21,500 1640 1640 871 871 522 

12 = 13 23 GIC/in.-lb/in.2

 GIIC/in.-lb/in.2 GIIIC/in.-lb/in.2 

0.30 0.45 0.486 3.17 3.17 1.8 

 

σ௡

଴/psi σ௦଴/psi σ௧

଴/psi K1/psi K2/psi K3/psi 

7800 12600 12600 1.37e9 1.37e9 1.37e9 

The middle strip of elements shown in blue in Figure 2 in the case of the FEAMAC simulations, are represented 

with a FEAMAC micromechanics repeating unit cell (RUC) through an ABAQUS user material. Usually such an 

RUC consists of a number of subcells and is used to represent a composite material. However, for these DCB 

simulations, the material is treated as homogenized and transversely isotropic, so the RUC consists of a single 

graphite/epoxy subcell with X-Z plane of isotropy, and interfacial debonding activated at the RUC interface as 

indicated in Figure 2. The same elastic transversely-isotropic graphite/epoxy material properties are used within the 

single-subcell FEAMAC RUC. It should be noted that FEAMAC enforces periodicity conditions at the RUC 

boundaries, thus the shown RUC represents a material that is continuous in the Y-direction, even though an 

interfacial displacement discontinuity, according to the ECI model described above, is enabled in the X-direction. 

The remainder of the geometry (yellow) is meshed using bilinear incompatible-mode plane strain (CPE4I) elements 

with elastic orthotropic material model with transverse isotropy as given in Table 1. 

In the case of simulations utilizing the ABAQUS cohesive and VCCT approaches, the geometry shown in 

Figure 2 remains the same with the exception of the introduction of a surface between the two plies. For simulations 

with the cohesive approach, cohesive surface elements are employed to discretize this surface between the two plies, 

whereas simulations utilizing the VCCT approach, this surface between the two plies is defined as the possible crack 

path. In both approaches, same crack length of CL is maintained and the remainder of the model is meshed using 

bilinear incompatible-mode plane strain (CPE4I) elements with the standard ABAQUS elastic orthotropic model 

with transverse isotropy, representing the same graphite/epoxy material discussed above.

B. Comparison of VCCT, Cohesive, and FEAMAC predictions for DCB Test 

A comparison of reaction forces predicted by the FEAMAC approach, VCCT approach, and cohesive surface 

approach is shown in Figure 3. In the case of VCCT simulation, debonding is determined entirely by the critical 

strain energy release rate GIC reported in Table 1. For cohesive surface approach, normal debond strength ௡

଴ as well 

L

CL

2W

Debond Fixed BC 

interface 

FEAMAC RUC 

Y 

X 

Applied displacements 0.05 in. 
NASA/TM—2010-216825 5 

Figure 3.—Reaction forces predicted by VCCT, cohesive, and 

FEAMAC using adjusted Gc and ௡

଴ for DCB. 

a) Cohesive Model b) VCCT Model c) FEAMAC Model 

Figure 4.—Comparison of stress component along the fiber (vertical) direction for cohesive, VCCT, 

and FEAMAC models at the beginning of crack propagation for a similar width of the geometry. 

as GIC have been used to determine the crack initiation and propagation, respectively. For the FEAMAC approach, 

this same debond strength ௡

଴ is used for the bond strength ஽஻

௡ parameter as expressed in Equation (2). In addition 

to the ஽஻

௡ , two other parameters critical to interfacial time-dependent behavior of the ECI model are the Λ௜ and ܤ௝

as reported in Equation (4). To create an interfacial debonding behavior similar to the cohesive and VCCT 

approaches, a value of 0.1 is selected for Λ௜ and 3.0 is selected for the parameter ܤ௝. 

The stress distributions for these three models along the crack front are shown in Figure 4. As expected with the 

denser cohesive mesh, the local stress distribution near the crack-tip is relatively high compared to the VCCT and 

FEAMAC DCB models. Both VCCT and FEAMAC models use the same element size of 0.01 in. over the entire 

length of the beam. Regardless of the disagreement of local peak stresses near the crack-tip, the overall stress 

distribution shows similarity while the overall response of these three models remain consistant as evidenced from 

reaction forces shown in Figure 3. 

FEAMAC Model 
NASA/TM—2010-216825 6 

TABLE 2.—TOTAL CPU TIME TAKE FOR COMPLETING 

A SIMULATION WITH FEAMAC, VCCT, 

AND COHESIVE MODELS 

Figure 5.—Geometry of the ENF specimen with applied boundary conditions (Ref. 15). 

The time taken to complete the execution of each model is shown in Table 2. Clearly, a significant execution 

time advantage is achieved when employing the FEAMAC approach over the VCCT and cohesive approaches. The 

DCB model with cohesive elements having a maximum allowable time of 0.01 sec, took approximately ten times 

more CPU time as compared to the FEAMAC simulation with a 0.075 sec constant time increment. If a larger 

maximum allowable time increment is used for the cohesive model, the solution process becomes unstable due to 

severe nonlinearities at the verge of the crack initiation. Abaqus/Standard offers a set of stabilization mechanisms to 

handle nonlinear problems by allowing the user to specify a number of solution control parameters such as number 

of equilibrium iterations (I0) and number of equilibrium iterations after which the logarithmic rate of convergence 

check begins (IR). Even with the current increment of 0.01 sec for the cohesive model, these convergence parameters 

had to be modified in order to stabilize the solution process. 

C. Simulations With an End-Notched Flexure Specimen 

The three-point end-notched flexure (ENF) test has been used to obtain the fracture toughness in Mode-II type 

delaminations where the initial crack is forced to grow by shearing. The dimensions and geometry used for ENF 

simulations is same as the DCB simulations, except for the applied boundary conditions (see Fig. 5), wherein one 

end is pinned and the other end is roller-supported. In contrast to the separation of laminates in the DCB test, a 

displacement at the middle of the specimen is used in the ENF simulations. 

Modeling of this geometry with FEAMAC, VCCT, and cohesive elements, is identical to the respective DCB 

simulations, that is the material data employed related to VCCT simulations (critical strain energy release rate GIIC

reported in Table 1 determines debonding) and for the cohesive surface approach, the tangential debond strength ௦

଴

as well as GIIC are used to determine the crack initiation and propagation, respectively. For the FEAMAC model, the 

same debond strength ௦

଴ has been used as the bond strength ஽஻

௧ as expressed in Equation (3). In addition to the 

஽஻

௧ , two other parameters critical to interfacial time-dependent behavior of the ECI model are the Λ௜ and ܤ ௝as 

reported in Equation (4). To create an interfacial debonding behavior similar to the cohesive and VCCT approaches, 

a value of 510–4 is selected for Λ௜ while 0.3 is selected for the parameter ܤ௝. 

D. Comparison of VCCT, Cohesive, and FEAMAC predictions for ENF Test 

To verify the material data presented in Table 1 related to Mode-II delamination, and identify the simulationspecific parameters associated with VCCT, FEAMAC, and cohesive surface approaches, ENF simulations have 

been conducted for comparison purposes. The reaction force measured at the center of the beam versus the applied 

displacement is shown in Figure 6. Clearly, all three models start debonding at approximately the same load level, 

yet the subsequent crack propagation pattern is distinctly different for the FEAMAC model as compared to both the 

VCCT and cohesive surface models. The latter two models qualitatively behave similarly with a sudden drop 

followed by a similar propagation pattern after the initiation of the crack. Note, the nature of damage progression in 

an ENF specimen is sudden and does not follow a slow crack propagation pattern as does the DCB simulation. Due 

to this sudden nature of crack propagation, it is essential to cognitive of the energies dissipated due to artificial 

damping introduced for numerical convergence purposes. Both cohesive and VCCT approaches needed significant 

Model Total CPU Time/s 

FEAMAC Model 1423 

VCCT Model 3241 

Cohesive Model 14,806 
NASA/TM—2010-216825 7 

Figure 6.—Reaction forces predicted by VCCT, cohesive, 

and FEAMAC using adjusted Gc and 0 n for ENF. 

Figure 7.—Initiation and propagation of a crack front between two layers. 

damping in order to maintain their numerical stability during simulations; whereas FEAMAC did not require any 

artificial damping or fine-tuning control parameters associated with the iterative convergence process. Therefore, in 

addition to execution time, the time needed to condition the solution process, i.e., experiment with numerous control 

parameters with FEAMAC is minimal, if not nonexistant, as compared to the VCCT and cohesive surface 

approches. 

E. Effect of Mesh Size on FEAMAC 

Finite element simulations of the damage and failure of composite materials have displayed a strong mesh 

dependency in crack initiation and propagation. This has been reflected on the global load-displacement response as 

well as the onset, size, and orientation of localized damage zones and crack growth in the composite materials. 

According to the implementation of FEAMAC, the stiffness of the entire element in a given direction reduces as the 

crack initiates at an integration point of an element, and progresses between two layers. Due to this stiffness 

degradation, a large crack opening (compared to a crack tip) is created and the mesh has to be fine enough to 

represent the propagating crack front (see Fig. 7). This mesh-dependence is common to other continuum damage 

mechanics approaches which rely upon stiffness degradation as their load redistribution mechanism. Consequently, 

to simulate the crack propagation as accurately as possible, the mesh must be refined along the possible crack front 

so that the debonding becomes more of a localized phenomenon, and the loss of element stiffness does not create an 

unrealistically large crack opening surface which would decrease the crack driving forces. The refined mesh used for 

FEAMAC is shown in Figure 8(a) together with the mesh used for VCCT simulations in which the mesh can be 

relatively coarser. 

 
NASA/TM—2010-216825 8 

 (a). FEAMAC Mesh (b). VCCT Mesh 

Figure 8.—Mesh refinement needed for FEAMAC compared to VCCT model. 

F. Uncoupled Normal and Shear in FEAMAC 

According to the ECI model in Equation (1), debonding can take place due to a normal or shear stress on any 

given face of a unit-cell. When one of these stresses exceeds the corresponding debonding strength, debonding 

initiates and the stiffness at an integration point has to degrade for both normal as well as shear. It was found that the 

current version of FEAMAC has not been implemented to account for this situation properly, in that if propogation 

of a crack were determined completely by the shear stress, the normal debond strength would still remain as shown 

in Figure 9(b), whereas conversely, if debonding initiates due to normal stress along the face of a unit-cell, it will 

still maintain a shear debond strength and the crack would continue to grow due to the normal debonding as in 

Figure 9(c). To overcome this situation a Hashin’s failure criterion has been used to couple the normal and shear 

stresses on a given face normal to the x2 or x3 direction as expressed in Equation (8): 

222 222

22 21 23 33 31 32

2 3 1 or 1 t t

nss nss

db db db db db db

F F

           

           

         

 (8) 

so that both normal and shear degradation of a given interface intiate simultaneously, if the stresses at the integration 

point satisfy the condition given in Equation (8). The simulations performed after this implementation of Hashin’s 

criteria displayed the desired effect of combined normal and shear failure as the crack continues to grow in a nonself-similar fashion. 

G. Prediction of Crack Path with FEAMAC 

In the previous sections, DCB and ENF simulations have been used to verify the material data for simulations 

under Mode-I and Mode-II loading conditions. Additionally, these simulations have been used to identify the 

appropriate simulation-specific conditions such as viscous regularization, mesh size, control parameters related to 

convergence, etc., for VCCT and cohesive simulations, as well as the Λ௜ and ܤ ௝parameters for the FEAMAC 

simulations. In this section, the ECI model available within FEAMAC, with the newly implemented Hashin’s 

criterion for interfacial failure, is used to demonstrate the prediction of crack path, when the crack path is unknown. 

To compare the FEAMAC predictions against other available techniques such as XFEM, simulation of a threepoint bend specimen under the plane strain conditions subjected to Mode-I loading, has been selected. The total 

length of the specimen is taken as 2 in. with a height of 0.4 in. (Ref. 17), and an initial crack length of 0.08 in. was 

employed as shown in Figure 10. 
NASA/TM—2010-216825 9 

(a) Eccentrically loaded 3-Point Bend Test 

(b) Crack propagation due to shear debonding (c) Crack propagation due to normal debonding 

Figure 9.—Shear and normal debonding within a 3-point bend specimen. 

Figure 10.—Three-point bend specimen with initial crack, boundary conditions, and loading. 

An angular rotation of 0.003 radians is applied at the center of the left and the right edges of the test specimen, 

while maintaining kinematic coupling constraints between each center point and the corresponding edge. Since the 

focus of this research is on the study of applying FEAMAC for predicting crack propagation with monolothic 

materials, isotropic material properties with damage parameters reported under Table 3 are used for these 

simulations. 

Simulation of the problem with VCCT, FEAMAC, and cohesive surface approaches are very much similar to 

the simulation of DCB under Mode-I loading using each technique. Unlike VCCT and cohesive approaches, XFEM 

and FEAMAC do not require a pre-existing crack or any prior knowledge of the crack path. Since the loading is 

symmetric on the three-point bending specimen, the crack is supposed to follow the expected crack path, as shown 

in the Figure 10. In the cohesive and VCCT models, this path has to be specified as the prospective crack path and 

the models will not look for other possible crack directions. However, in the XFEM simulation, crack direction can 

change based on the specified crack propagation criteria. Any purturbation in the numerical solution process may 

cause the crack to change its direction, even though the loading, boundary conditions, geometry, and meshing are 

perfectly symmetrical (Ref. 17). Therefore, to compare the results with the VCCT and cohesive approaches, it is 

necessary to maintain the symmetry of displacement by constraining the displacements of the left and right edges 

where the angular displacement is applied. This precaution is necessary for XFEM simulations to force the 

propagation of the crack in the anticipated direction. In fact, the need to utilize constraint equations to maintain 

symmetry is a deficiency in the implementation of XFEM, since the crack should have propagated along the path 

highlighted in Figure 10 without introducing additional constraints. 

Angular Rotation 

0.003 rad 

Angular Rotation 

0.003 rad 

Bonded Surface (Expected Crack Path) 

Initial Crack 
NASA/TM—2010-216825 10 

TABLE 3.—MATERIAL PROPERTIES USED FOR THREE-POINT BEND SIMULATIONS 

E/ksi  GIC/in.-lb/in.2

 GIIC/in.-lb/in.2 GIIIC/in.-lb/in.2 

29,000 0.3 0.571 0.571 0.571 2.284 

௡

଴ /ksi ௦଴/ksi ௧

଴/ksi 

25.4 25.4 25.4 

Only for FEAMAC Only for XFEM 

ௗ௕

௡ /ksi ௗ௕௦ /ksi Λ௜ ܤ ௝Maximum Principal Stress/ksi 

33.4 33.4 0.10 0.30 18.9 

Figure 11.—Reaction moment predicted by FEAMAC, 

VCCT, Cohesive, and XFEM simulations. 

The reaction moment observed at the center point of the right-edge verses the applied angular displacement is 

plotted in Figure 11 for all considered techniques. Reaction values predicted at the initiation of the crack are similar 

among all the techniques, but the reaction values observed during the propagation of the crack are somewhat 

different among the models. Both VCCT and XFEM have sudden crack opening, while cohesive and FEAMAC 

predict relatively slower crack extension against the applied angular rotation. 

As expected, the crack in the VCCT and cohesive models followed the predefined crack path along the center of 

the specimen (indicated as the bonded surface in Figure 10). Since the symmetric conditions have been enforced 

during XFEM simulation, the crack continues to grow along the middle of the specimen, similar to the VCCT and 

cohesive simulations. It is worth mentioning that the addition of such constraint equations to an FEA model can 

significantly increase the computational time, especially with nonlinear material properties and large-scale structural 

problems. The displacement discontinuity approach in FEAMAC does not require additional constraints to maintain 

symmetry, and it was observed that the crack continued along the center similar to the three above techniques. Based 

on the results reported in Figure 11, it can be concluded that the ability to predict the crack path, and the resultant 

reaction moment from each technique, have been verified through the simulations with three-point bend specimen. 

 
NASA/TM—2010-216825 11 

Figure 12.—Eccentrically loaded three-point bend specimen with initial crack, boundary conditions, and loading. 

(a) Crack path and plot of failed with XFEM (b) Crack path and failed elements with 

FEAMAC together with crack predicted 

by XFEM (in red)

Figure 13.—Comparison of crack path predicted by XFEM and FEAMAC. 

The finite element-based VCCT or cohesive element approaches require that the path of the crack be known in 

advance, whereas FEAMAC and XFEM have the capability to predict where the crack is going to grow without any 

prior knowledge about the crack path. To demonstrate this capability with FEAMAC, an eccentrically loaded 

3-point bend specimen with the same dimensions and initial crack length as in the case described above, is used with 

displacement and boundary conditions as shown in Figure 12. 

Both the XFEM and FEAMAC approaches are used to obtain the crack propagation through finite element 

simulations. Predicted crack patterns for the given displacement are shown in Figure 13(a) for XFEM and 

Figure 13(b) for FEAMAC. Results shown in Figure 13(a) display the XFEM predicted crack path (shown in white) 

together with the partially failed elements as the crack propagates through the beam, whereas Figure 13(b) shows the 

the XFEM crack path (in red) superimposed on the FEAMAC elements that failed completely during the crack 

propagation. These two figures demonstrate a fundamental difference between the XFEM and FEAMAC 

approaches, in that the XFEM approach allows an actual crack surface to be generated by the crack front whereas 

theFEAMAC crack path is described by the elements that become near zero-stiffness elements to mimic the 

propagating crack, creating a property discontinuity. Therefore, FEAMAC will have the typical mesh sensitivity 

exhibited by continuum damage mechanics approaches where no physical discontinuites are present. 

Displacement of 0.0016 in.

Initial Crack length of 0.08 in. 
NASA/TM—2010-216825 12 

According to these simulation results, both approaches predict the crack path closely at the beginning of the 

crack propagation, but they deviate as the crack progresses along the specimen. Additionally, FEAMAC tends to fail 

a number of elements across the width of the crack path whereas XFEM fails only the elements that are directly 

divided by the propagating crack front. This discrepency could be associated with the debonding parameters utlized, 

the use of the Hashin failure criteria, and/or mesh sensitivity issues; all of which will be the subject of future work. 

Selection of an appropriate principal stress value for XFEM simulations has a significant effect on the prediction of 

the crack path as well as the extent of the the crack for a given displacement. Adequate mesh refinement is another 

factor to obtain a reliable prediction for crack path with XFEM. Since the entire element fails as the crack 

progresses, FEAMAC also relies on a finer mesh density to maintain a narrow crack path through the material. 

IV. Conclusions 

An extensive comparison of the displacement discontinuity approach in FEAMAC with existing approaches 

(cohesive elements, the Virtual Crack Closure Technique, and XFEM method) available in Abaqus finite element 

software has been conducted to evaluate the micromechanics-based approach and identify needed enhancements. As 

a result of these simulations, issues related to the implementation of FEAMAC as well as conceptual issues such as 

shear and normal coupling at debonding, have been identified and resolved. During this study, the FEAMAC 

approach has been used to investigate failure within materials efficiently and to explore the possibility of predicting 

non-self-similar crack paths. An eccentrically loaded 3-point bend test has been simulated with FEAMAC and the 

results are similar to XFEM predictions at the beginning of crack initiation, with some deviation from the XFEM 

predictions as the crack front progresses through the specimen under displacement control. 

Based on the experience with finite element based approaches and the micromechanics approach with 

FEAMAC, it is understood that the micromechanics based FEAMAC approach has far less sensitivity to parameters 

associated with numerical convergence. Additionally, FEAMAC offers significant advantage over the number of 

elements used compared to the cohesive approach and hence the computational time required to complete a 

simulation. This may offer an advantage over established finite element-based approaches in providing efficient and 

robust design tools for bonded composite joints. Most significantly, FEAMAC appears to have the potential to 

predict arbitrary crack paths similar to XFEM, whereas the VCCT and cohesive element approaches require 

predetermined crack paths. 

To validate fully the FEAMAC simulation results, it is necessary to compare the predictions with experimental 

data. In the absence of such data, comparison has been made to an XFEM simulation of an isotropic, eccentricallyloaded, three-point bend specimen. Isotropic material properties where utilized due to the difficulty in estimating a 

suitable principal stress value for XFEM, which prevented using orthotropic material properties. In contrast 

FEAMAC simulations need only the debonding strength values that have been verified through DCB and ENF 

simulations and thus easily applied to orthotropic and composite materials. 

References 

1. Bednarcyk, B.A. and Arnold, S.M., “A New Local Debonding Model with Application to the Transverse 

Tensile and Creep Behavior of Continuously Reinforced Titanium Composites,” NASA/TM—2000-210029, 

2000. 

2. Bednarcyk, B.A. and Arnold, S.M., “Transverse Tensile and Creep Modeling of Continuously Reinforced 

Titanium Composites with Local Debonding,” International Journal of Solids and Structures 39 (7), 2002, 

pp. 1987–2017. 

3. Bednarcyk, B.A., and Arnold, S.M., MAC/GMC 4.0 User’s Manual, NASA/TM—2002-212077, 2002. 

4. Krueger, R., “The Virtual Crack Closure Technique: History, Approach and Applications,” NASA/CR—2002-

211628, 2002. 

5. Comanho, P.P. and Davila, C.G., “Mixed-Mode Decohesion Finite Elements for the Simulation of 

Delamination in Composite Materials,” NASA/TM—2002-211737, 2002. 

6. Ranatunga, V., Bednarcyk, B.A. and Arnold, S.M., “Modeling Delamination in composites via Continuum 

Interfacial Displacement Discontinuities,” 50th AIAA ASME/ASCE/AHS/ASC SDM Conference, Palm 

Springs, CA, May, 2009. 

7. Moes, N., Dolbow, J., Belytschko, T., “A Finite Element Method for Crack Growth Without Remeshing,” 

International Journal for Numerical Methods in Engineering, 46(1), 1999, pp. 131–150. 

8. Jones, J.P., and Whittier, J.S., “Waves at Flexibly Bonded Interfaces,” Journal of Applied Mechanics 34, 1967, 

pp. 905–909. 
NASA/TM—2010-216825 13 

9. Aboudi, J., “Damage in Composites – Modeling of Imperfect Bonding,” Composites Science and Technology 

28, 1987, pp. 103–128. 

10. Sankurathri, A., Baxter, S. and Pindera, M.J., “The Effect of Fiber Architecture on the Inelastic Response of 

Metal Matrix Composites with Interfacial and Fiber Damage,” in Damage and Interfacial Debonding in 

Composites, G.Z. Voyiadjis and D.H. Allen (eds.), Elsevier, New York, 1996, pp. 235–257. 

11. Achenbach, J.D. and Zhu, H., “Effect of Interfacial Zone on Mechanical Behavior and Failure of FiberReinforced Composites,” Journal of the Mechanics and Physics of Solids, 37 (3), 1989, pp. 381–393. 

12. Wilt, T.E. and Arnold, S.M., “Micromechanics Analysis Code (MAC) User Guide: Version 2.0,” NASA TM 

107290, 1996. 

13. Needleman, A., “A Continuum Model for Void Nucleation by Inclusion Debonding,” Journal of Applied 

Mechanics 54, 1987, pp. 525–531. 

14. Robertson, D.D. and Mall, S., “Micromechanical Analysis of Metal Matrix Composite Laminates with 

Fiber/Matrix Interfacial Debonding,” Composites Engineering 4 (12), 1994, pp. 1257–1274. 

15. Song, K., Davila, C.G., and Rose, C.A., “Guidelines and Parameter Selection for the Simulation of Progressive 

Delamination,” Abaqus Users’ Conference, Newport, RI, 2008. 

16. Rudraraju S., Salvi, A., Garikpati, K., and Waas, A., “In-Plane fracture of laminated fiber reinforced composites 

with varying fracture resistance: experimental observations and numerical crack propagation simulations,” 50th 

AIAA ASME/ASCE/AHS/ASC SDM Conference, Palm Springs, CA, May, 2009. 

17. Dassault Systemes, “Modeling Fracture and Failure with Abaqus,” Lowell, MA, 2009. 
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-10-2010 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Modeling Progressive Damage Using Local Displacement Discontinuities Within the 

FEAMAC Multiscale Modeling Framework 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Ranatunga, Vipul; Bednarcyk, Brett, A.; Arnold, Steven, M. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 645846.02.07.03.13 

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

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-2010-216825 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Category: 39 

Available electronically at http://gltrs.grc.nasa.gov 

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

13. SUPPLEMENTARY NOTES

14. ABSTRACT

A method for performing progressive damage modeling in composite materials and structures based on continuum level interfacial 

displacement discontinuities is presented. The proposed method enables the exponential evolution of the interfacial compliance, resulting in 

unloading of the tractions at the interface after delamination or failure occurs. In this paper, the proposed continuum displacement 

discontinuity model has been used to simulate failure within both isotropic and orthotropic materials efficiently and to explore the possibility 

of predicting the crack path, therein. Simulation results obtained from Mode-I and Mode-II fracture compare the proposed approach with the 

cohesive element approach and Virtual Crack Closure Techniques (VCCT) available within the ABAQUS (ABAQUS, Inc.) finite element

software. Furthermore, an eccentrically loaded 3-point bend test has been simulated with the displacement discontinuity model, and the 

resulting crack path prediction has been compared with a prediction based on the extended finite element model (XFEM) approach.

15. SUBJECT TERMS

Advanced composites; Structural systems; High temperature; FEAMAC 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

19 

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