Robert K. Goldberg and Kelly S. Carney

Glenn Research Center, Cleveland, Ohio

Paul DuBois

George Mason University, Fairfax, Virginia 

Canio Hoffarth and Subramaniam Rajan

Arizona State University, Tempe, Arizona

Gunther Blankenhorn

Livermore Software Technology Corporation, Livermore, California

Incorporation of Plasticity and Damage Into an 

Orthotropic Three-Dimensional Model With 

Tabulated Input Suitable For Use in 

Composite Impact Problems

NASA/TM—2015-218849

August 2015
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated 

to the advancement of aeronautics and space science. 

The NASA Scientific and Technical Information (STI) 

Program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Officer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI Program provides access 

to the NASA Technical Report Server—Registered 

(NTRS Reg) and NASA Technical Report Server—

Public (NTRS) thus providing one of the largest 

collections of aeronautical and space science STI in 

the world. Results are published in both non-NASA 

channels and by NASA in the NASA STI Report 

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of significant 

scientific and technical data and information 

deemed to be of continuing reference value. 

NASA counter-part of peer-reviewed formal 

professional papers, but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

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

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis.

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA.

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 757-864-6500

• Telephone the NASA STI Information Desk at

757-864-9658

• Write to:

NASA STI Program

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199
Robert K. Goldberg and Kelly S. Carney

Glenn Research Center, Cleveland, Ohio

Paul DuBois

George Mason University, Fairfax, Virginia 

Canio Hoffarth and Subramaniam Rajan

Arizona State University, Tempe, Arizona

Gunther Blankenhorn

Livermore Software Technology Corporation, Livermore, California

Incorporation of Plasticity and Damage Into an 

Orthotropic Three-Dimensional Model With 

Tabulated Input Suitable For Use in 

Composite Impact Problems

NASA/TM—2015-218849

August 2015

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

30th Annual Technical Conference

sponsored by the American Society for Composites

East Lansing, Michigan, September 28–30, 2015
Acknowledgments

Authors Hoffarth and Rajan gratefully acknowledge the support of the Federal Aviation Administration through 

Grant Number 12-G-001 entitled “Composite Material Model for Impact Analysis,” William Emmerling, Technical Monitor.

Available from

Trade names and trademarks are used in this report for identification 

only. Their usage does not constitute an official endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

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

NASA STI Program

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22161

703-605-6000

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

Incorporation of Plasticity and Damage Into an Orthotropic 

Three-Dimensional Model With Tabulated Input Suitable 

For Use in Composite Impact Problems

Robert K. Goldberg and Kelly S. Carney

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Paul DuBois

George Mason University

Fairfax, Virginia 22030

Canio Hoffarth and Subramaniam Rajan

Arizona State University

Tempe, Arizona 85281

Gunther Blankenhorn

Livermore Software Technology Corporation

Livermore, California 94551

Abstract

The need for accurate material models to simulate the deformation, damage and failure of polymer 

matrix composites under impact conditions is becoming critical as these materials are gaining increased 

usage in the aerospace and automotive industries. While there are several composite material models 

currently available within commercial transient dynamic finite element codes, several features have been 

identified as being lacking in the currently available material models that could substantially enhance the 

predictive capability of the impact simulations. A specific desired feature pertains to the incorporation of 

both plasticity and damage within the material model. Another desired feature relates to using

experimentally based tabulated stress-strain input to define the evolution of plasticity and damage as 

opposed to specifying discrete input properties (such as modulus and strength) and employing analytical 

functions to track the response of the material. To begin to address these needs, a combined plasticity and 

damage model suitable for use with both solid and shell elements is being developed for implementation 

within the commercial code LS-DYNA. The plasticity model is based on extending the Tsai-Wu 

composite failure model into a strain-hardening based orthotropic plasticity model with a non-associative 

flow rule. The evolution of the yield surface is determined based on tabulated stress-strain curves in the 

various normal and shear directions and is tracked using the effective plastic strain. The effective plastic 

strain is computed by using the non-associative flow rule in combination with appropriate numerical 

methods. To compute the evolution of damage, a strain equivalent semi-coupled formulation is used, in 

which a load in one direction results in a stiffness reduction in multiple coordinate directions. A specific 

laminated composite is examined to demonstrate the process of characterizing and analyzing the response 

of a composite using the developed model.
NASA/TM—2015-218849 2

Introduction

As composite materials are gaining increased use in aircraft components where impact resistance 

under high energy impact conditions is important (such as the turbine engine fan case), the need for 

accurate material models to simulate the deformation, damage and failure response of polymer matrix 

composites under impact conditions is becoming more critical. Within commercially available transient 

dynamic finite element codes such as LS-DYNA (Ref. 1), there are several material models currently 

available for application to the analysis of composites. The available models include relatively simple

models such as a derivative of the Chang-Chang model (Ref. 2), where criteria related to ratios of stresses 

to failure strengths are used to signify failure, and the composite elastic constants are selectively reduced 

based on the failure mode. More sophisticated material models incorporated within transient dynamic 

finite element codes include continuum damage mechanics based models such as the model developed by 

Matzenmiller, et al. (Ref. 3). In this approach, the initiation and accumulation of damage is assumed to be 

the primary driver of any nonlinearity in the composite response. The failure stresses and strains of the 

material in each of the coordinate directions are specified by the user, the evolution of the damage is 

computed based on “damage variables”, and the nonlinearity of the material stress-strain response is 

approximated based on the input data and the evolution of the damage. While not necessarily included 

within commercial codes currently, another approach that has been used to model the strain rate 

dependent response of a composite is to assume that all of the nonlinearity is due to deformation 

mechanisms. An example of this approach was developed by Sun and Chen (Ref. 4), where a quadratic 

plastic potential function was developed and the plastic strains were computed based on the gradient of 

the plastic potential function. Stress-strain curves obtained for various fiber orientation angles are used to 

characterize the coefficients in the plastic potential function based on the values required to collapse the 

various separate stress-strain curves into a master curve.

While the material models discussed above and other models have been utilized with some level of 

success in modeling the nonlinear and impact response of polymer composites, there are some areas 

where the predictive capability can be improved. Most importantly, the existing models often require 

significant a priori knowledge of the damage and failure in the analyzed structure such that their use as 

predictive tools can be limited. While these models generally assume that the nonlinear response of the 

composite is due to either deformation mechanisms (such as plasticity) or damage mechanisms, an 

improved model should have the capability to simulate the actual material behavior in which the material 

nonlinearity is due to a combination of both deformation and damage mechanisms. The input to current 

material models generally consists of point-wise properties (such as a specified failure stress or failure 

strain) that lead to curve fit approximations to the material stress-strain curves. This type of approach 

leads either to models with only a few parameters, which provide a crude approximation at best to the 

actual stress-strain curve, or to models with many parameters which require a large number of complex 

tests to characterize. An improved approach would be to use tabulated data from a well-defined set of 

experiments to accurately define the complete stress-strain response of the material. Furthermore, many 

of the existing models are only suitable for use with two-dimensional shell elements, which cannot 

capture the through-thickness response, which may be significant in impact applications. Ideally, a fully 

three-dimensional formulation suitable for use with solid elements would be desirable, along with a shell 

element formulation.

To begin to address these needs, a multi-institution consortium has been formed to develop and 

implement a new composite material model within LS-DYNA, which will be implemented as MAT_213.

The material model is meant to be a fully generalized model suitable for use with any composite 

architecture (unidirectional, laminated or textile). For the deformation model, the commonly used Tsai-
NASA/TM—2015-218849 3

Wu composite failure criteria has been generalized and extended to a strain-hardening model with a 

quadratic yield function and a non-associative flow rule. The coefficients of the yield function for the new 

composite model are determined based on tabulated stress-strain curves in the various normal and shear 

directions, along with selected off-axis curves. The non-associative flow rule is used to compute the 

components of the plastic strain along with the effective plastic strain. The evolution of the yield stresses 

in the various coordinate directions is tracked based on the current value of the effective plastic strain. For 

the damage model, a strain equivalent formulation has been developed, which allows the plasticity and 

damage calculations to be uncoupled, and the plasticity calculations to take place in the effective 

(undamaged) stress space. In traditional damage mechanics models such as the one developed by 

Matzenmiller et al. (Ref. 3), a load in a particular coordinate direction is assumed to result in a stiffness 

reduction only in the direction of the applied load. However, for reasons to be discussed later in this 

paper, in the current model a semi-coupled formulation is developed in which a load in one direction 

results in a stiffness reduction in all of the coordinate directions.

In the following sections of this paper, a summary of the derivation of the rate-independent 

deformation model is presented. The procedures to be used to characterize the material constants in the 

yield function and the flow law are discussed. Next, the rationale for, and detailed derivation of, the semicoupled damage model is discussed along with a summary of the test matrix that will be required to 

properly characterize and validate the developed model. Finally, selected verification studies that have 

been conducted to ensure that the deformation model has been implemented correctly are presented.

Deformation Model Derivation

A general quadratic three-dimensional orthotropic yield function based on the Tsai-Wu failure model 

is specified as follows, where 1, 2, and 3 refer to the principal material directions.

( ) ( )

( )









































σ

σ

σ

σ

σ

σ









































+ σ σ σ σ σ σ









































σ

σ

σ

σ

σ

σ

σ = +

31

23

12

33

22

11

66

55

44

13 23 33

12 22 23

11 12 13

11 22 33 12 23 31

31

23

12

33

22

11

1 2 3

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0

0 0 0

0 0 0

0 0 0

F

F

F

F F F

F F F

F F F

f a F F F

(1)

In the yield function, σij represents the stresses and Fij and Fk are coefficients that vary based on the 

current values of the yield stresses in the various coordinate directions. By allowing the coefficients to 

vary, the yield surface evolution and hardening in each of the material directions can be precisely defined.

The values of the normal and shear coefficients can be determined by simplifying the yield function for 

the case of unidirectional tensile and compressive loading in each of the coordinate directions along with 

shear tests in each of the shear directions, with results as shown below
NASA/TM—2015-218849 4

2

31

66

33 33

33

33 33

3

2

23

55

22 22

22

22 22

2

2

12

44

11 11

11

11 11

1

1 1 1 1

1 1 1 1

1 1 1 1

1

σ = σ σ = σ − σ =

σ = σ σ = σ − σ =

σ = σ σ = σ − σ =

= −

F F F

F F F

F F F

a

T C T C

T C T C

T C T C

(2)

In the above equation, the stresses are the current value of the yield stresses in the normal and shear 

directions (determined using procedures to be discussed below), where the superscript T indicates the 

tensile yield stress, and the superscript C indicates the absolute value of the compressive yield stress. To 

determine the values of the off-axis coefficients (which are required to capture the stress interaction 

effects), the results from 45° off-axis tests in the various coordinate directions can be used. An important 

point to note is that due to experimental or numerical variability, or alternatively just due to the 

fundamental behavior of the material, computing the off-diagonal terms of the yield function in this 

manner may result in a yield function that is not convex (which is a requirement for plasticity theory 

(Ref. 5)). As a result, to satisfy the requirements of the chosen yield function, the off-diagonal terms may 

need to be adjusted based on the values of the other coefficients in the yield function in order to ensure 

convexity of the yield surface.

A non-associative flow rule is used to compute the evolution of the components of plastic strain. The 

plastic potential for the flow rule is shown below

2 66 312 55 232 12 11 22 23 22 33 31 33 11 44 122 33 332 22 222 11 11 σ + σ + σ + 2 σ σ + 2 σ σ + 2 σ σ + σ + σ + σ

=

H H H H H H H H H

h

(3)

where σij are the current values of the stresses and Hij are independent coefficients, which are assumed to 

remain constant. The values of the coefficients are computed based on average plastic Poisson’s ratios 

(Ref. 6). The plastic potential function in Equation (3) is used in a flow law to compute the components 

of the plastic strain rate, where the usual normality hypothesis from classical plasticity (Ref. 5) is assumed 

to apply and the variable λ

 is a scalar plastic multiplier.

σ

ε

∂

∂ = λ h

p   (4)

Given the flow law, the principal of the equivalence of plastic work (Ref. 5) can be used to determine 

expressions for the effective stress and effective plastic strain. By following this procedure, one can 

conclude that the plastic potential function h can be defined as the effective stress and the plastic 

multiplier can be defined as the effective plastic strain rate.

To compute the current value of the yield stresses needed for the yield function, the common practice 

in plasticity constitutive equations is to use analytical functions to define the evolution of the stresses as a 

function of the components of plastic strain (or the effective plastic strain). Alternatively, in the 

developed model tabulated stress-strain curves are used to track the yield stress evolution. The user is 

required to input twelve stress versus plastic strain curves. Specifically, the required curves include 

uniaxial tension curves in each of the normal directions (1,2,3), uniaxial compression curves in each of 

the normal directions (1,2,3), shear stress-strain curves in each of the shear directions (1-2, 2-3 and 3-1), 
NASA/TM—2015-218849 5

and 45° off-axis tension curves in each of the 1-2, 2-3 and 3-1 planes. The 45° curves are required in 

order to properly capture the stress interaction effects. By utilizing tabulated stress-strain curves to track 

the evolution of the deformation response, the experimental stress-strain response of the material can be 

captured exactly without any curve fit approximations. The required stress-strain data can be obtained 

either from actual experimental test results or by appropriate numerical experiments utilizing stand-alone 

codes. Currently, only static test data is considered. Future efforts will involve adding strain rate and 

temperature dependent effects to the computations. To track the evolution of the deformation response 

along each of the stress-strain curves, the effective plastic strain is chosen to be the tracking parameter.

Using a numerical procedure based on the radial return method (Ref. 5) in combination with an iterative 

approach, the effective plastic strain is computed for each time/load step. The stresses for each of the 

tabulated input curves corresponding to the current value of the effective plastic strain are then used to 

compute the yield function coefficients.

Damage Model Derivation

The deformation portion of the material model provides the majority of the capability of the model to 

simulate the nonlinear stress-strain response of the composite. However, in order to capture the nonlinear 

unloading and local softening of the stress-strain response often observed in composites (Ref. 7), a 

complementary damage law is required. In the damage law formulation, strain equivalence is assumed, in 

which for every time step the total, elastic and plastic strains in the actual and effective (undamaged) 

stress spaces are the same (Ref. 8). The utilization of strain equivalence permits the plasticity and damage 

calculations to be uncoupled, as all of the plasticity computations can take place in the effective stress 

space.

The first step in the development of the damage model is to relate the actual stresses to a set of 

effective stresses by use of a damage tensor M

σ = Mσeff (5)

The effective stress rate tensor can be related to the total and plastic strain rate tensors by use of the 

standard elasto-plastic constitutive equation

( ) σeff C ε ε p  =  −  (6)

where C is the standard elastic stiffness matrix and the actual total and plastic strain rate tensors are used 

due to the strain equivalence assumption. By differentiating Equation (5) and substituting in Equation (6), 

the actual stress rate can be written in terms of the total and plastic strain rates and the actual stress. In the 

following expression, Voigt notation is assumed to be appropriate and the damage tensor is assumed to be 

invertible.

σ MC(ε ε ) MM σ

σ Mσ Mσ

−1 = − +

= +

   

  

p

eff eff (7)

An algorithm to carry out the uncoupled plasticity/damage analysis is summarized below. In the 

algorithm, the superscript “n” represents values computed in the previous time step, and the superscript 

“n+1” indicates values to be computed in the current time step. In the first step of the algorithm, the actual 

stresses are converted into effective stresses using the damage tensor M. In the second step, the plasticity 
NASA/TM—2015-218849 6

calculations are carried out in the effective stress space to compute the current value of the plastic strain 

rate, and the effective stress values are updated. Next, in step 3 the damage tensor is modified based on 

the computed plastic strain rate. Finally, in step 4 the modified damage tensor is used to compute the 

updated values of the actual stresses based on the updated effective stresses. The algorithm is summarized 

symbolically below, where Δt is the time step.

( )

( )

[ ]

1 1 1

1

1

1

4.

3.

2.

1.

+ + +

+

+

−

=

= + ∆

= + − ∆

=

n

eff n n

p n n

p n

eff neff

n n n

eff

t

σ M σ

M M M ε

σ σ C ε ε

σ M σ



  (8)

To ensure that the strain equivalence assumption is valid for this approach, the equations utilized in the 

algorithm should be able to be manipulated to yield an equation similar to that shown in the second 

expression in Equation (7), which was derived based solely on the strain equivalence assumption. By 

starting with the equation shown in step 4 of the algorithm, substituting the equation used in step 2 of the 

algorithm in for the modified effective stress, and breaking up the modified damage tensor into the sum of 

the current value of the damage tensor and the increment in the damage tensor (ΔM), the following set of 

equations is developed.

( ) ( )

( ) ( )

( ) ( ) ( ) ( )

( ) ( )

( ) ( ) n n

p n n n

n n n p n n

n n n p n n

n n n p n n

p n n n n

t t

t

t

t

t

M σ

M M C ε εσ σ

σ M C ε ε σ M M σ

σ M M C ε ε M M M σ

σ M C ε ε M M σ

σ M M σ C ε ε

1 1

1

1 1 1

1 1

1 1 1 1

1 1 1

+ − +

+ + −

+ −

+ + + −

+ + −

∆

∆ = − +

∆

−

= − ∆ + + ∆

= + ∆ − ∆ + + ∆

= − ∆ +



  = + − ∆

 

 

 

 

 

(9)

To a first order approximation, the last expression in Equation (9) corresponds to the second expression in 

Equation (7), which indicates that strain equivalence assumption is appropriate for the combined 

plasticity/damage algorithm.

Given the usual assumption that the actual stress tensor and the effective stress tensor are symmetric, 

the actual stresses can be related to the effective stresses in the following manner, where the damage 

tensor M is assumed to have a maximum of 36 independent components.

[ ]









































σ

σ

σ

σ

σ

σ

=









































σ

σ

σ

σ

σ

σ

eff

xz

eff

yz

eff

xy

eff

zz

eff

yy

eff

xx

xz

yz

xy

zz

yy

xx

M (10)
NASA/TM—2015-218849 7

In many damage mechanics models for composites, for example (Refs. 3 and 7), the damage tensor is 

assumed to be diagonal or manipulated to be a diagonal tensor, leading to the following form.

[ ]









































=

66

55

44

33

22

11

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

M

M

M

M

M

M

M (11)

The shear terms M44, M55 and M66 can be independent (such as in Matzenmiller, et al. (Ref. 3)) or 

functions of the normal damage terms M11, M22 and M33.

The implication of a diagonal damage tensor is that loading the composite in a particular coordinate

direction only leads to a stiffness reduction in the direction of the load due to the formation of matrix 

cracks perpendicular to the direction of the load. However, several recent experimental studies 

(Refs. 9, 10, and 11) have shown that in actual composites, particularly those with complex fiber 

architectures, a load in one coordinate direction can lead to stiffness reductions in multiple coordinate 

directions. One example of this phenomena can be seen in Figure 1, which came from research conducted 

by Salem and Wilmoth (Ref. 11). The figure shows a triaxially braided composite which has been loaded 

in tension in the transverse direction and then reloaded in compression in the longitudinal direction. As 

can be seen in the figure, after loading in the transverse direction significant damage above and beyond 

simple matrix cracks was present in the material. As a result, when the material was reloaded in the 

longitudinal direction, the measured longitudinal modulus was significantly reduced from the baseline 

value, indicating that the damage resulting from the transverse load affected the stiffness in the 

longitudinal direction. There have been limited attempts to incorporate this damage coupling into an 

analytical technique, for example by Voyiadjis and Park (Ref. 12) and Bednarcyk, et al. (Ref. 13).

However, these efforts were developed within the overall concept of a damage mechanics theory in which 

all of the nonlinearity of the composite response was assumed to be due to damage mechanisms.

Furthermore, in the developed theoretical approaches analytical functions were used to track the evolution 

of the damage and the reduction of the stiffness. The efforts described in the current paper are geared 

towards developing a damage theory, uncoupled from the plasticity theory, which tracks the stiffness 

reduction and damage accumulation as a function of the plastic strain by the use of tabulated input.

One approach to incorporating the coupling of damage modes would be to use a non-diagonal damage 

tensor, such as the one shown below for the case of plane stress.

























σ

σ

σ





















=





















σ

σ

σ

eff

xy

eff

yy

eff

xx

xy

yy

xx

M M M

M M M

M M M

61 62 66

21 22 26

11 12 16

(12)

However, while this formulation would allow for directional coupling, it would have the side effect of a 

unidirectional load in the actual stress space resulting in a multiaxial load in the effective stress space. For 

the strain equivalent combined plasticity-damage formulation envisioned for this model, this would be an 

undesirable side effect as the plasticity calculations could be adversely affected due to the introduction of 

nonphysical stresses.
NASA/TM—2015-218849 8

Figure 1.—Damage observed in triaxially braided composite 

subjected to transverse loading followed by longitudinal 

loading.

Figure 2.—Coupled damage resulting from loading in x coordinate direction.

To avoid the undesired stress coupling, a diagonal damage tensor is required. However, to account for 

the damage interaction in at least a semi-coupled sense, each term in the diagonal damage matrix should 

be a function of the plastic strains in each of the normal and shear coordinate directions, as follows for the 

example of the M11 term for the plane stress case

( ) p

xy

p

yy

p M11 = M11 εxx ,ε ,ε (13)

To explain this concept graphically, a schematic is shown in Figure 2 for the case of loading in the x 

coordinate direction. A plastic strain is applied to an undamaged specimen, with an original area Axx

perpendicular to the x axis and an original area Ayy perpendicular to the y axis. The undamaged modulus 

in the x direction is Exx and the undamaged modulus in the y direction is equal to Eyy. The specimen is 

damaged due to the plastic strain. The original specimen is unloaded and reloaded elastically in the x 
NASA/TM—2015-218849 9

direction. Due to the damage, the reloaded specimen has a reduced area in the x direction of dxx Axx and a 

reduced modulus in the x direction of dxx Exx . The reduced area and modulus are a function of the damage 

induced by the plastic strain in the x direction as follows

( ( )) ( ( )) xx p

xx xxxx dxxxx

xx pxx xxxx dxxxx

A d A

E d E

= − ε

= − ε

1

1

(14)

where xx dxx is the damage in the x direction due to a load in the x direction. Alternatively, if the damaged 

specimen was reloaded elastically in the y direction, due to the assumed damage coupling the reloaded 

specimen would have a reduced area in the y direction of dxx Ayy and a reduced modulus in the y direction 

of dxx Eyy due to the load in the x direction. The reduced area and modulus are again a function of the 

damage induced by the plastic strain in the x direction as follows

( ( ))

( ( )) yy

p

xx yyxx dxxyy

yy

p

xx yyxx dxxyy

A d A

E d E

= − ε

= − ε

1

1

(15)

where yy dxx is the damage in the y direction due to a load in the x direction. Similar arguments can be 

made and equations developed for the situation where the original specimen is loaded plastically in the y 

direction.

The next issue is how to properly model the damage coupling for the case of multiaxial loading. A 

simple case for consideration is a composite specimen being simultaneously strained plastically in the x 

and y directions. One way to approach the problem would be to assume that the loss in area in the x 

direction due to the straining in the x direction is as follows.

xx xx xx dxx Axx − Axx = A d (16)

Next, one can assume that the loss in area in the x direction due to the loading in the y direction is as 

follows

xx xx yy

dyy Axx − Axx = A d (17)

where dyy Axx is the reduced area in the x direction resulting from a load in the y direction and xx dyy is the 

damage in the x direction resulting from a load in the y direction. The total loss in area in the x direction 

could be computed by adding together the two separate area losses to obtain a total loss in area.

[ ] xx

yy xx xx xx dyyxx dxx 2Axx − Axx − A = A d + d (18)

The ratio of the total damaged area to the undamaged area in the x direction can now be computed, which 

theoretically would lead to the total amount of damage in the x direction in the composite under biaxial 

loading
NASA/TM—2015-218849 10

[ ] xx

yy xxxx

xx

xx

yy xx xx xx xx

xx

dyy

xx dxx xx xx xx d d

A

A A d d

A

A A A A = − − − + = − + + 1 2 (19)

The relationship between the actual stress and effective stress in the x direction would then be the 

following.

( ( ) ( )) eff

xx p

yy xxyy

p

xx xx σxx = 1− dxx ε − d ε σ (20)

The error in this approach is that the assumption is made that the load in the x and y directions are 

both acting on undamaged areas. In reality, the loads are acting on damaged areas, and just adding to the 

damaged area. For example, if one loaded the material in the y direction first, the reduced area in the x 

direction would be equal to dyy Axx and the reduced modulus in the x direction would be equal to dyy Exx . If 

one would then subsequently load the material in the x direction, the baseline area in the x direction 

would not equal the original area Axx, but the reduced area dyy Axx . Likewise, the baseline modulus in the x 

direction would not be equal to the original modulus Exx, but instead the reduced modulus dyy Exx .

Therefore, the loading in the x direction would result in the following further reduction in the area and 

modulus in the x direction.

( ( )) ( ( )) ( ( ))

( ( )) ( ( )) ( ( )) xx p

yy xxyy

p

xx xxxx dyyxx pxx xxxx dxxxx

xx pyy xxyy

p

xx xxxx dyyxx pxx xxxx dxxxx

A d A d d A

E d E d d E

= − ε = − ε − ε

= − ε = − ε − ε

1 1 1

1 1 1

(21)

These results suggest that the relation between the actual stress and the effective stress should be based on 

a multiplicative combination of the damage terms as opposed to an additive combination of the damage 

terms. For example, for the case of plane stress, the relation between the actual and effective stresses 

could be expressed as follows

( )( )( )

( )( )( )

( )( )( ) eff

xy

xy

xy

xy

yy

xy

xy xx

eff

yy

yy

xy

yy

yy

yy

yy xx

eff

xx xxxy xxyy xx xx xx

d d d

d d d

d d d

σ = − − − σ

σ = − − − σ

σ = − − − σ

1 1 1

1 1 1

1 1 1

(22)

where for each of the damage terms the subscript indicates the direction of the load which initiates the 

particular increment of damage and the superscript indicates the direction in which the damage takes 

place. Note that for the full three-dimensional case the stress in a particular coordinate direction is a 

function of the damage due to loading in all of the coordinate directions (x, y, z, xy, xz and yz). By using 

a polynomial to describe the damage, the coupled terms represent the reduction to the degree of damage 

resulting from the fact that in a multiaxial loading case the area reductions are combined.

To properly characterize the damage model, an extensive set of test data is required. Due to the 

tabulated nature of the input, each of the damage parameters ( xx dxx , yy dxx , etc.) has to be determined as a 

function of the plastic strain in a particular coordinate direction (such as p εxx ). For example, to determine 

the damage terms for the case of loading in the x direction, a composite specimen has to be loaded to a 

certain plastic strain level in the x direction. The material is then unloaded to a state of zero stress, and 

then reloaded elastically in each of the coordinate directions to determine the reduced modulus of the 
NASA/TM—2015-218849 11

material in each of the coordinate directions. Expressions such as those in Equations (14) and (15) can 

then be used to determine the required damage parameters for the particular value of plastic strain. The

process needs to be repeated for multiple values of plastic strain in the x direction in order to establish the 

full characterization of the variation of the damage parameters as a function of the plastic strain in the x 

direction.

Verification Studies For Deformation Model

A set of verification studies for the deformation portion of the material model were conducted using 

data for a T800S/3900-2B unidirectional composite (Ref. 14). This particular composite is composed of 

intermediate modulus, high strength fibers embedded within a toughened epoxy matrix. Full details of the 

verification studies will be provided in a future paper, however a summary is provided here for 

completeness. An important point to note is that at the current time only the deformation portion of the 

material model has been implemented numerically within the LS-DYNA computer code, so only 

deformation analyses will be discussed here. Future efforts will involve implementing and verifying the 

damage theory described above. However, since the evolution of the parameters in the damage model will 

be based on the plastic strains, ensuring that the nonlinear deformation response of the composite is being 

properly simulated by the material model is critical to ensuring that the subsequent characterization and 

verification of the damage model is correct.

The input data utilized for the verification studies is a combination of actual experimental data obtained 

by Raju and Acosta (Ref. 14) and numerical simulations. The numerical simulations were conducted to 

obtain required stress-strain curves which were not available from the provided experimental data. To 

conduct the numerical simulations, micromechanics analyses were conducted where the properties of the 

fiber and matrix were used to simulate the overall response of the composite. The micromechanics analyses 

were conducted by using a combination of purely analytical simulations conducted using the NASA Glenn 

developed MAC/GMC code based on the Generalized Method of Cells (Ref. 15) and finite element models 

where the fiber and matrix were explicitly simulated. The fiber was assumed be linear elastic with 

transversely isotropic properties and the matrix was assumed to be isotropic with an elastic-plastic material 

response. A fiber volume fraction of 0.54 was assumed based on data presented in Bogert, et al (Ref. 16). 

The constituent material properties were chosen to correlate with composite experimental data obtained by 

Raju and Acosta (Ref. 14) and Bogert, et al. (Ref. 16). The elastic modulus and Poisson’s ratio of the matrix 

were chosen to be values representative of epoxy based materials (Ref. 17). The fiber longitudinal modulus, 

transverse modulus and longitudinal Poisson’s ratio were chosen such that the composite longitudinal 

modulus, transverse modulus and longitudinal Poisson’s ratio computed by the micromechanics analyses 

correlated with the values obtained by Bogert, et al. (Ref. 16). For the fiber transverse Poisson’s ratio, a 

value representative of that utilized for similar carbon fibers (Ref. 17) was used. The fiber in-plane shear 

modulus was chosen such that the computed composite in-plane shear modulus correlated to the composite 

in-plane shear modulus determined by Raju and Acosta (Ref. 14). The transverse shear moduli were 

computed by appropriate assumptions based on the assumed transverse isotropy of the unidirectional 

material. The yield stress for the matrix was selected such that the computed in-plane shear stress-strain 

curve for the composite roughly correlated with the experimental stress-strain curve determined by Raju and 

Acosta (Ref. 14). Full details of the numerical experiments and the process used to correlate the properties 

for the numerical experiments will be provided in a future paper. The correlated fiber and matrix properties 

used for the analyses are provided in Table 1.

To conduct the verification studies, finite element models such as the ones shown in Figure 3 for 

tension and Figure 4 for shear were used. Note that 64 eight noded solid elements were used for the 
NASA/TM—2015-218849 12

analyses. For the tension simulations, the nodes on the left hand side of the model were constrained and a 

constant displacement in the x direction was applied to the nodes on the right face. For the shear 

simulations, the nodes on the bottom surface were constrained and a displacement in the x direction was 

applied to the top surface.

TABLE 1.—FIBER AND MATRIX CONSTITUTIVE 

PROPERTIES USED FOR MICROMECHANICS ANALYSES

Property Fiber Matrix

E11,GPa 275.6 3.45

E22, GPa 15.5 15.5

ν12 0.20 0.35

ν23 0.25 0.35

G12, GPa 103.4 1.27

σyield, MPa N/A 137.8

Figure 3.—Finite element model for tension verification analyses.

Figure 4.—Finite element model for shear verification analyses.
NASA/TM—2015-218849 13

Simulated stress-strain curves were computed for a variety of load cases, examples of which are 

shown in Figure 5 for the case of [0°] tension, Figure 6 for the case of in-plane shear in the x-y plane and 

Figure 7 for [45°] off-axis tension in the x-y plane. In all of the figures the “Experimental” curves were 

the curves that were provided as input data for the material model which were based on the actual or 

numerically generated stress-strain curves, and the “Simulated” curves were the curves computed using 

the material model. In all cases, the simulated curves correlated very well to the input curves, indicating 

that the material model was able to accurately represent the nonlinear deformation of the composite in 

terms of correctly replicating the input data. Future efforts will involve conducting more complex 

validation analyses which will include analyzing the deformation response of more complex laminated 

composites as well as simulating the deformation response of a composite under impact conditions.

Figure 5.—Verification analyses for 0° tension test.
NASA/TM—2015-218849 14

Figure 6.—Verification analyses for in-plane shear test.

Figure 7.—Verification Analyses for 45° off-axis tension test.
NASA/TM—2015-218849 15

Conclusions

A generalized composite model suitable for use in polymer composite impact simulations has been 

developed. The theory for the rate independent deformation and damage portions of the composite model 

have been developed and numerical implementation of the deformation model has been completed. The 

complete composite model will be implemented into the LS-DYNA commercial transient dynamic finite 

code as MAT 213. For the deformation model, the Tsai-Wu composite failure model has been generalized 

into an orthotropic yield function with a non-associative flow rule. Tabulated stress-strain data is utilized 

to track the deformation response of the material, using the effective plastic strain as the tracking variable.

A strain equivalent damage model has been developed in which loading the material in a particular 

coordinate direction can lead to damage in multiple coordinate directions. The actual and effective 

stresses are related by multiplicative combinations of the various damage variables.

Future efforts will involve generalizing the deformation model to incorporate the ability to simulate 

the effects of strain rate and temperature on the material response. The damage model will be numerically 

implemented within the LS-DYNA code. Methods to model failure and element removal will be 

developed and implemented into LS-DYNA. Extensive additional verification and validation studies will 

be conducted to verify the accuracy and capability of the overall material model. Overall, when 

completed the material model when implemented into MAT 213 will provide significant improvements to 

the state of the art in the modeling of the impact response of polymer composites.

References 

1. Hallquist, J. 2013. LS-DYNA Keyword User’s Manual, Version 970. Livermore Software Technology 

Corporation, Livermore, CA.

2. Chang, F.-K. and K.-Y. Chang. 1987. “A Progressive Damage Model for Laminated Composites 

Containing Stress Concentrations,” Journal of Composite Materials, 21:834-855.

3. Matzenmiller, A., J. Lubliner, and R.L. Taylor. 1995. “A constitutive model for anisotropic damage in 

fiber-composites,” Mechanics of Materials, 20:125-152.

4. Sun, C.T., and J.L. Chen. 1989. “A Simple Flow Rule for Characterizing Nonlinear Behavior of Fiber 

Composites,” Journal of Composite Materials, 23:1009-1020.

5. Khan, A.S., and S. Huang. 1995. Continuum Theory of Plasticity. John Wiley and Sons, New York.

6. Goldberg, R., K. Carney, P. DuBois, C. Hoffarth, J. Harrington, S. Rajan, and G. Blankenhorn. 2014. 

“Theoretical Development of an Orthotropic Elasto-Plastic Generalized Composite Model,” 

NASA/TM—2014-218347, National Aeronautics and Space Administration, Glenn Research Center, 

Cleveland, OH.

7. Barbero, E.J. 2013. Finite Element Analysis of Composite Materials Using ABAQUS. CRC Press, 

Boca Raton, FL.

8. Lemaitre, J., and R. Desmorat. 2005. Engineering Damage Mechanics: Ductile, Creep and Brittle 

Failures. Springer, Berlin.

9. Ogasawara, O, T. Ishikawa, T. Yokozeki, T. Shiraishi, and N. Watanabe. 2005. “Effect of on-axis 

tensile loading on shear properties of an orthogonal 3D woven SiC/SiC composite,” Comp. Sci. 

Technol., 65:2541-2549.

10. Salavatian, M., and L.V. Smith. 2014. “The effect of transverse damage on the shear response of fiber 

reinforced laminates,” Comp. Sci. Technol., 95:44-49.

11. Salem, J., and N. Wilmoth, 2015. Personal Communication.

12. Voyiadjis, G.Z., and T. Park. 1995. “Anisotropic Damage of Fiber-Reinforced MMC Using Overall 

Damage Analysis,” Journal of Engineering Mechanics, 121(11):1209-1217.

13. Bednarcyk, B.A., B. Stier, J.-W. Simon, S. Reese, and E.J. Pineda. 2015. “Meso- and micro-scale 

modeling of damage in plain weave composites,” Composite Structures, 121:258-270.
NASA/TM—2015-218849 16

14. Raju, K.S., and J.F. Acosta. 2010. “Crashworthiness of Composite Fuselage Structures—Material 

Dynamic Properties, Phase I,” DOT/FAA/AR-09/8, U.S. Department of Transportation, Federal 

Aviation Administration, Washington, D.C.

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

NASA/TM—2002-212077/VOL2, National Aeronautics and Space Administration, Glenn Research 

Center, Cleveland, OH.

16. Bogert, P.B., A. Satyanarayana, and P.B. Chunchu. 2006. “Comparison of Damage Path Predictions 

for Composite Laminates by Explicit and Standard Finite Element Analysis Tools,” 47th 

AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference,

American Institute for Aeronautics and Astronautics, Washington, D.C.

17. Murthy, P.L.N., C.A. Ginty, and J.G. Sanfeliz, 1993. “Second Generation Integrated Composite 

Analyzer (ICAN) Computer Code,” NASA TP–3290, National Aeronautics and Space 

Administration, Lewis Research Center, Cleveland, OH.