Robert Goldberg and Kelly Carney

Glenn Research Center, Cleveland, Ohio

Paul DuBois

George Mason University, Fairfax, Virginia

Canio Hoffarth, Joseph Harrington, and Subramaniam Rajan

Arizona State University, Tempe, Arizona

Gunther Blankenhorn

Livermore Software Technology Corporation, Livermore, California 

Theoretical Development of an Orthotropic ElastoPlastic Generalized Composite Material Model

NASA/TM—2014-218347

August 2014
	
	
	








	



	

!

	
	
"#"	




		


	
	

		


"$
	

%	



	&	

		


'	
	
	



	
(

	)
	

			


	




	
*
	

)	

+

")
		*



*	
	",

. 5$6$9:;9$%)
	




		
	
<
	 


				


	
	
=
	
	

"






	







		


		
	+	*
	


	

		





	=


		


. 5$6$9>5>%)(:> 

		
	"
	 


&	?#	 

	
	*
	#	
	






(




="

. $%)$%))5%)

"+

	 


	
		

. $%@5)5$5:;9$%$


		




	"


	
	
	




	
	


	"

. 5$9:;9$%


	
	
	

	

 

	
	
	
<





	*<

	

. 5$6$9)9%5+

	


	


		
'



&	
	



	

&
	&

			

@
	

	
	




	
	


*,

. 	
	




http://www.sti.nasa.gov

. 5+
"
	?

help@sti.nasa.gov

. @="
	?




	

(#BBFHIJIHJKLF

. 

	

(#

 BBFHIJIHJKLQ

. R	
,


	

(#

 $	
		

	



IUUJ	(	

6
	>(QULIVHUFQL
Robert Goldberg and Kelly Carney

Glenn Research Center, Cleveland, Ohio

Paul DuBois

George Mason University, Fairfax, Virginia

Canio Hoffarth, Joseph Harrington, and Subramaniam Rajan

Arizona State University, Tempe, Arizona

Gunther Blankenhorn

Livermore Software Technology Corporation, Livermore, California 

Theoretical Development of an Orthotropic ElastoPlastic Generalized Composite Material Model

NASA/TM—2014-218347

August 2014


	



	


W)	$	 

$%
BBUFJ

		
	

UF:	
9+(X:	$
	



	"9	

	
*	

"$
	
	


(	
	>YKHULQLUB
Acknowledgments


	6
	6		
)<	"#
*
	
@	

	
	


W	ZUQ+W+LLU[$


>	>

	
"\R
5

	>

	

	



$	
		

	



IUUJ	(	

6
	>(QULIVHUFQL



	

	

JFLU*)


=	]QQFUQ

	
",^^***


	
	
	#		
	
	



"	




	


	=	
	
"
	



	


Level of Review,
	"	*"


NASA/TM—2014-218347 1 

Theoretical Development of an Orthotropic Elasto-Plastic 

Generalized Composite Material Model 

Robert Goldberg and Kelly Carney 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Paul DuBois 

George Mason University 

Fairfax, Virginia 22030 

Canio Hoffarth, Joseph Harrington, and Subramaniam Rajan 

Arizona State University 

Tempe, Arizona 85004 

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 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 LSDYNA (Livermore Software Technology Corporation), there are several features that have been 

identified that could improve the predictive capability of a composite model. To address these needs, a 

combined plasticity and damage model suitable for use with both solid and shell elements is being 

developed and is being implemented into LS-DYNA as MAT_213. A key feature of the improved 

material model is the use of tabulated stress-strain data in a variety of coordinate directions to fully define 

the stress-strain response of the material. To date, the model development efforts have focused on 

creating the plasticity portion of the model. The Tsai-Wu composite failure model has been generalized 

and extended to a strain-hardening based orthotropic yield function with a nonassociative flow rule. The 

coefficients of the yield function, and the stresses to be used in both the yield function and the flow rule, 

are computed based on the input stress-strain curves using the effective plastic strain as the tracking 

variable. The coefficients in the flow rule are computed based on the obtained stress-strain data. The 

developed material model is suitable for implementation within LS-DYNA for use in analyzing the 

nonlinear response of polymer composites. 

Introduction 

As composite materials are gaining increasing 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 LS-DYNA (Ref. 1), several 

material models are available for application to the analysis of composites. For example, the ChangChang failure model (Ref. 2) is utilized in MAT_22 and MAT_54. In these models, combinations of 

ratios of stresses to failure strengths are utilized to predict fiber or matrix based failure. The response is 
NASA/TM—2014-218347 2 

assumed be linear elastic, with limited capability to simulate the nonlinear shear response. In MAT_22 

the failure is assumed to be brittle, while in MAT_54 the composite elastic constants are selectively 

reduced based on the failure mode, and a gradual unloading is permitted until ultimate element failure is 

reached. In MAT_58, a continuum damage model developed by Matzenmiller et al (Ref. 3) is employed, 

where the initiation and accumulation of damage is assumed to be the primary driver of nonlinearity in 

the composite response. The failure stresses and strains of the material in each of the coordinate directions 

are specified by the user, and the material stress-strain curves are approximated based on this data. The 

original version of the material model assumes a material response that is independent of strain rate. 

However, in the enhanced model MAT_158, a modified version of the model used in MAT_58 is 

employed which incorporates strain rate dependence into the material response (Ref. 1). A viscoelastic 

Prony series based on shear moduli is used to modify the computed stresses. The strain rate dependence in 

all of the coordinate directions is assumed to be identical. MAT_161 is an additional continuum damage 

mechanics model available within LS-DYNA (Ref. 4). In this model, a stress to strength ratio approach in 

a variety of coordinate directions is used to specify the beginning of either fiber or matrix based failure. 

Appropriate failure mode based damage functions are then used to compute the reduction in elastic 

moduli in each of the coordinate directions. In MAT_219, based on the CODAM model developed at the 

University of British Columbia (Ref. 5), a strain versus failure strain ratio based approach is used to 

predict the initiation of damage in a sublaminate of the composite. Separate damage accumulation and 

modulus reduction functions based on failure mode and coordinate direction are implemented in a 

continuum damage mechanics formulation. In MAT_221, damage accumulation functions based on 

current strains, damage initiation strains and failure strains are used to reduce the elastic moduli of the 

composite in each of the coordinate directions in a damage mechanics approach (Ref. 1). In MAT_261, 

which is based on a damage and failure model developed by Pinho et al (Refs. 6 and 7), separate models 

for fiber tension failure, fiber kinking failure, matrix tensile failure and matrix compressive failure are 

developed, where different functional forms are used for each of the failure criteria (for example, the 

matrix compression failure criterion is based on an extension of the Mohr-Coulomb failure criterion). The 

various failure criteria are combined together using fracture mechanics concepts to develop a constitutive 

model. In MAT_262, based on a model developed by Camonho, et al (Refs. 8 and 9), an energy approach 

is utilized to generate damage functions in various coordinate directions within the context of a 

continuum damage mechanics formulation. 

While all of the material models discussed above have been utilized with some level of success in 

modeling the impact response of polymer composites, there are some areas where the predictive 

capability can be improved. In general, the existing models often require significant a priori knowledge of 

the damage and failure 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 (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 the current material models generally consists of point-wise properties that lead 

to curve fit approximations to the material stress-strain curves. This type of approach leads to either 

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 throughthickness response which may be significant in impact applications. Ideally, a fully three-dimensional 

formulation suitable for use with solid elements would be available, along with the shell element 

formulation. In addition, incorporating the effects of strain rate is only possible in a few of the existing 
NASA/TM—2014-218347 3 

models, and in those models, the effects of strain rate on the material response are assumed to be the same 

in all of the material coordinate directions. 

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. 

Currently, the primary focus of the effort has been on the development and testing of the orthotropic 

macroscopic three-dimensional plasticity model. As will be discussed in more detail in the following 

sections of this paper, the commonly used Tsai-Wu composite failure criterion (Ref. 10) has been 

generalized and extended to a strain-hardening model with a quadratic yield function and a nonassociative 

flow rule. A similar approach has been used by Sun and Chen (Ref. 11), in which a general quadratic plastic 

potential function was generalized into a plasticity model suitable for use with composites. 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 stresses that 

are used in the flow rule are computed based on the input tabulated stress-strain curves. Systematic 

procedures have been developed to compute the various coefficients in the yield function and the flow rule 

based on the tabulated input data. An important point to note is that the developed material law is not limited 

to the analysis of unidirectional plies alone. Alternatively, the material model is meant to be a fully 

generalized model suitable for use with any composite architecture (unidirectional, laminated or textile). In 

the following sections of this paper, a detailed derivation of the yield function and the flow law are 

presented. The procedures to be used to characterize the material constants in the yield function and the flow 

law are also discussed in detail. Details of the numerical implementation of the material model within LSDYNA and descriptions of verification and validation studies that were conducted using the developed 

material model are available in a companion paper (Ref. 12). 

Material Law 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. 
NASA/TM—2014-218347 4 

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. 

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. For example, 

the stresses in the local material axis system resulting from a uniaxial tensile test of a [45°] composite (for 

a unidirectional composite), or any uniaxial tensile test conducted at 45° in the 1-2 plane from the 

longitudinal (1) material axis for a multiply laminated or textile composite, can be computed to be the 

following by the stress transformation equations (Ref. 10). 

12 45

22 45

11 45

0.5

0.5

0.5

Y

Y

Y


  
 



  



  


 (3) 

In the above equation, Y45 is the structural level current uniaxial yield stress resulting from a tensile test 

on a [45°] ply, and the remaining stresses are the material axis system based stresses. By substituting this 

equation into Equation (1) and solving for the off-axis constant F12, the following expression can be 

obtained 

  11 22 44

45

1 2

2

45

12

2

2 1 F F F F F F

Y Y


  




 
 
  (4) 

where all of the terms are as defined earlier. By using similar procedures involving 45° rotations in the 2-

3 plane and the 1-3 plane, similar expressions for the constants F23 and F13 can be determined. 

A nonassociative 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

31 66 223 55 2

12 12 11 22 23 22 33 31 33 11 44 2

33 33 2

22 22 2

11 11
  
  
  2 
 
  2 
 
  2 
 
  
  
  




H H H H H H H H H

h

 (5) 

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

remain constant. Procedures for determining the values of the coefficients will be discussed later in this 

paper. The plastic potential function in Equation (5) is used as follows in a flow law, where the usual 
NASA/TM—2014-218347 5 

normality hypothesis from classical plasticity (Ref. 13) is assumed to apply and the variable 

is a scalar 

plastic multiplier. 







66 31 12

55 23 23

44 12 12

13 11 23 22 33 33 33

12 11 22 22 23 33 22

11 11 12 22 13 33 11

2

2

2

2 2 2

2

2 2 2

2

2 2 2

2






 






 






 


  
  




 


  
  




 


  
  




 



  

H

h

H

h

H

h

H H H

h

H H H

h

H H H

h

h

p

p

p

p

p

p















p

 (6) 

In the above equation, p

ij are the components of plastic strain rate. For the shear components, tensorial, 

not engineering, strains are used. By examining ratios of the components of plastic strain in response to a 

given uniaxial stress component, the normal and off-axis coefficients Hij in the plastic potential function 

can be related to “plastic Poisson’s ratios”, p

ij  , as follows. This concept, as will be shown later in this 

paper, is useful in developing procedures to characterize the values of the coefficients. 

33

13

33

11

31

22

23

22

33

23

11

13

11

33

13

33

23

33

22

32

22

12

22

11

21

11

12

11

22

12

11 0 22 0 33 0

H

H

H

H

H

H

H

H

H

H

H

H

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

 
 

  
   
 

  
   
 



  


 
 

  
   
   
   
 

  



  
  
 























(7) 

In the above equation the expressions in the first column are only valid for the case where the only 

nonzero stress is in the 1 direction, the expressions in the second column are only valid for the case where 

the only nonzero stress is in the 2 direction, and the expressions in the third column are only valid for the 

case where the only nonzero stress in in the 3 direction. An important point to note is that while the yield 

function can accommodate differences between the tensile and compressive responses of the composite, 

the flow law, due to its lack of linear terms, cannot make such a distinction. However, introduction of a 

linear term would make the plastic Poisson’s ratio, and thus the coefficients in the flow law, highly 

dependent on the stress which would most likely lead to erratic behavior of the model. Another 

consideration is that the coefficients of the flow law could be assumed to vary based on the current stress 

and strain state, but that would require developing evolution conditions on the flow law, which would be 

very difficult to determine based on experimental data. 

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

determine expressions for the effective stress and effective plastic strain. By taking the vector product of 
NASA/TM—2014-218347 6 

the stress tensor and the plastic strain rate tensor, one ends up with the plastic potential function 

multiplied by the plastic multiplier 

. Since, by the principal of the equivalence of plastic work, the 

vector product of the stress tensor and the plastic strain rate tensor must equal the product of the effective 

stress and the effective plastic strain rate, 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. This concept is expressed mathematically in the equation shown below 

p e e p p h h W    
 



    



:  : (8) 

where e is the effective stress and p eis 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), 

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. For certain composite architectures, the number of stress-strain curves that need to be input 

can be simplified. For example, for unidirectional composites, the normal response in the 3 direction can 

be assumed to be identical to the normal response in the 2 direction due to the transverse isotropy of the 

material. Likewise, the shear response in the 1-3 direction can be assumed to be identical to the shear 

response in the 1-2 direction. Similarly, the 45° off-axis curve in the 1-3 plane can be assumed to be 

identical to the 45° off-axis curve in the 1-2 plane. Due to the transverse isotropy of a unidirectional 

composite, one can also demonstrate that for the case of a unidirectional composite the 45° off-axis curve 

in the 2-3 plane is approximately equal to the normal stress-strain curve in the 2 (or 3) direction. For other 

composite architectures with significant symmetry (such as two-dimensional plain or satin weave 

composites), a certain degree of simplification in the stress-strain data that is required to be input is also 

possible. 

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 the numerical procedure described briefly later 

in this paper and in more detail in a companion paper (Ref. 12), 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. To enable this 

tracking, each of the input stress-strain curves needs to be converted into a plot of stress versus effective 

plastic strain. To enable this conversion, the principal of the equivalence of plastic work is used again to 

relate the unidirectional plastic strain increment to the effective plastic strain increment. For example, for 

the case of unidirectional loading in the 1 direction the effective plastic strain increment (and integrated 

current total effective plastic strain) can be computed as follows 
NASA/TM—2014-218347 7 

 



 

 





 

 


   
 


   


   
   

11 11

11 11

11 11

11 11

11 11 11 11

11 11

H

d

H

d

d

d H d

d hd

d d hd

p

p

e

p

p

e

p

e p

p

e p

p

e p e e p

ij ij

 (9) 

where 11 is the unidirectional stress in the loading direction, p d 11  is the plastic strain increment in the 

loading direction and p de is the increment in effective plastic strain. Note that in a general case the 

increment in effective plastic strain can be computed as the incremental area under the stress-plastic strain 

curve divided by the current value of the effective stress. A graphical representation of this process is 

shown below in Figure 1. 

To compute the evolution of the effective plastic strain, a numerical algorithm based on the radial 

return method is employed. Given the revised value of the effective plastic strain, the yield stress values 

and the overall stress state of the material can be updated. Details of the numerical implementation are 

provided in a companion paper (Ref. 12), but some of the key theoretical fundamentals of the algorithm 

are given below. The standard elastic constitutive equation is used to compute the revised stresses, where 

the flow law is applied for the computation of the plastic strains. 

 



 



	



  
  
 



C   C 

h p : : (10) 

In the above equation, C is the standard elastic stiffness matrix, is the total strain rate tensor, and the 

remaining terms are as defined previously. To compute the effective plastic strain rate, the consistency 

condition is applied in combination with the elastic constitutive equation 

Figure 1.—Conversion of stress versus plastic strain curves to stress versus effective 

plastic strain curves. 


11

11

11 11

11

p

E




   



11

User provided load curves are

true stress versus true plastic strain

Internally stored is true stress

versus effective plastic strain

  h

p d e

11 11 
 


NASA/TM—2014-218347 8 

 

















 

 



  



 



	



 
 



 

 







 

d

f h f d

f

d

f h f d f

f f f

q

q

C



C 



q

q

C  C



q

q





:

:

: : 0

0







 (11) 

where f is the yield function defined in Equation (1) and q is the vector of yield stresses based on the 

input stress-strain curves. 

 12 23 31 45 45 23 45 31 11 22 33 11 22 33  
 
 
 
 
 
 
 
 
 
Y 
Y 
 
Y 
 T T T C C C q (12) 

where 
Y 45
23 is the yield stress resulting from a uniaxial tensile test conducted at 45° in the 2-3 plane 

and 
Y 45
31 is the yield stress resulting from a uniaxial tensile test conducted at 45° in the 3-1 plane. To 

compute the derivative of the components of the q vector with respect to the effective plastic strain , a 

chain rule type of computation is carried out where the derivatives of the individual components of the 

vector (the respective yield stresses) with respect to the effective plastic strain are determined as follows. 








  


d

d

d

d

d

d p

ij

p

ij

ij ij (13) 

The derivative of the components of the yield stress vector with respect to plastic strain is the 

instantaneous slope of the stress versus plastic strain curve from each of the 12 input stress-strain curves. 

The derivative of the plastic strain with respect to the plastic multiplier  can be computed from the flow 

law for each of the specialized unidirectional (or off-axis) cases represented by the input stress-strain 

curves. Note that for the off-axis cases the flow law must first be written in the material axis system and 

then the plastic strains must be transformed into the structural axis system in order to compute the proper 

derivative of the plastic strain with respect to the plastic multiplier. A secant iteration approach is used in 

combination with the radial return method to compute the specific value of the effective plastic strain rate 

(Ref. 12). 

Characterization of Flow Law Coefficients 

The coefficients in the flow rule need to be characterized based on data obtained from experimental 

stress-strain curves for the composite. Specifically, data from the 12 experiments leading to the input 

stress-strain curves must be utilized. The data can be obtained either from actual experiments or 

numerical experiments. For example, if one knows the mechanical properties and stress-strain curves for 

the fiber and matrix constituents, numerical analyses can be conducted using either high fidelity finite 

element analysis (where the fiber and matrix are modeled as distinct constituents) or analytical tools such 

as the NASA Glenn developed micromechanics code MAC/GMC (Ref. 14) to obtain numerical 

approximations to the actual composite experimental data. 

First, the procedures used to obtain the required coefficients for the simplified case of a unidirectional 

carbon fiber based polymer matrix composite will be discussed. Afterward, a generalized procedure to 

obtain the required constants for a general composite will be presented. For a unidirectional carbon fiber 

composite, as discussed by Sun and Chen (Ref. 11) a reasonable approximation to make is that due to the 
NASA/TM—2014-218347 9 

linear behavior of the carbon fiber the plastic strain in the fiber direction (the 1 direction) is equal to zero 

for all values of stress. By examining the second expression in Equation (6), one can conclude that for a 

general stress state the only way the plastic strain in the 1 direction can be zero is for the coefficients H11, 

H12 and H13 to all be identically zero. Next, since the transverse tension response of the composite can 

display a certain degree of nonlinearity, a reasonable approximation to make is that for the case of a 

unidirectional load in the 2 direction the effective stress is equal to the applied transverse tension stress 

22. By simplifying the plastic potential function h (shown in Equation (5)) for the case of a unidirectional 

transverse tension stress, and by setting this term (along with the transverse tensile stress) equal to the 

effective stress, one can show that the value of the coefficient H22 must equal one. Since a unidirectional 

composite is transversely isotropic, one can set the value of H33 to also equal one. Given these values, 

using the definition of the plastic Poisson’s ratio p

23  shown in Equation (7), one can find that the 

coefficient H23 equals the negative of the plastic Poisson’s ratio p

23  . Since the flow surface coefficients 

are assumed to be constant, an average value of the plastic Poisson’s ratio must be determined from the 

experimental and/or numerical test data from a unidirectional transverse tension test. To determine the inplane shear coefficient H44, utilizing a procedure similar to that developed by Sun and Chen (Ref. 11) for 

their model, the plastic potential function h given in Equation (5) and the plastic strain definition given in 

Equation (6) can be simplified for the case of in-plane shear loading in the 1-2 direction, yielding the 

following expressions. 

44

12

44 12

12 44

0.5

0.5

H

H

h H

p

p

e pe p

e



     

 
  


 (14) 

Starting from the experimental stress versus plastic strain plot for the case of in-plane shear loading, 

using the expressions in Equation (14) the value of H44 can be optimized to provide the closest possible 

match with the overall effective stress versus effective plastic strain curve, selected for the case of a 

unidirectional carbon fiber composite to be the transverse tension stress versus plastic strain curve. Due to 

the transverse isotropy of the unidirectional composite, the value of the coefficient H66 can be set equal to 

the value of the coefficient H44. To determine the value of the H55 coefficient, a similar procedure can be 

carried out for the case of shear loading in the 2-3 direction. Alternatively, due to the transverse isotropy 

of the composite, by computing the effective stress for the case of off-axis loading in the 2-3 plane, an 

isotropic type of relation can be determined for the coefficient H55

 p H55 23  2 1  (15) 

where p

23  is the plastic Poisson’s ratio in the 2-3 direction, obtained from the material data. 

Similar procedures are used for the case of a general composite with no identified symmetries. An 

example of a general composite would be a triaxially braided composite which is modeled as a smeared, 

homogenous material. Another example would be a generalized laminated composite modeled as a 

smeared, homogeneous material. For the case of a general composite, a unidirectional tension test in the 

longitudinal (1) direction is chosen as the baseline case. By utilizing Equation (5) and setting the plastic 

potential function and the unidirectional longitudinal stress equal to the effective stress, the value of the 

H11 coefficient can be found to equal one. Based on the relations given in Equation (7), the value of the 

coefficient H12 can be found to equal the negative of the plastic Poisson’s ratio p

12  , and the value of the 
NASA/TM—2014-218347 10 

coefficient H13 can be found to be equal to the negative of the plastic Poisson’s ratio p

13  . As before, the 

value of the plastic Poisson’s ratio used in the computations is the average value of the parameter 

determined from the experimental and/or numerically obtained data. In a similar fashion, the values of the 

coefficients H22, H23 and H33 can be computed. 

p

p

p

p

p

p

p

H

H

H

31

13 33

21

12

23 23

21

12 22



 



  




 

 (16) 

To compute the values of the coefficients H44, H55, and H66, the same procedure that was utilized to 

compute H44 for the case of a unidirectional carbon fiber composite is used, only all three shear loading 

1223, 31) must be used. 

Conclusions 

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

developed. The plasticity model will be part of a generalized plasticity and damage model which will be 

implemented into LS-DYNA as MAT_213. The yield function is based on the Tsai-Wu composite failure 

model, and a suitable nonassociative flow rule was defined. Methods of utilizing tabulated stress-strain 

data to track the evolution of the yield stresses as a function of the effective plastic strain have been 

developed. The elastic constitutive equation can be used in combination with the consistency condition to 

compute the effective plastic strain. In a companion paper (Ref. 12), further details of the numerical 

algorithm are provided along with examples of verification and validation studies which have been 

performed to examine the accuracy and efficiency of the analytical model. 

Future efforts will involve the development of an equivalent damage model, to be employed in 

combination with the plasticity model, which will also be based on utilizing tabulated stress-strain data to 

track the evolution of damage in the composite. Element failure will also be incorporated within the 

model. In addition, strain rate and temperature effects will be added to the analysis capability. Overall, 

when completed the composite model as implemented in MAT_213 will provide significant 

improvements to the state of the art in the modeling of the impact response of polymer matrix composites. 

References 

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

Corporation, Livermore, CA, 2013. 

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

Containing Stress Concentrations,” Journal of Composite Materials, Vol. 21, pp. 834-855, 1987. 

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

fiber-composites,” Mechanics of Materials, Vol. 20, pp. 125-152, 1995. 

4. Yen, C.F.: “Ballistic Impact Modeling of Composite Materials,” 7th International LS-DYNA Users 

Conference, Livermore Software Technology Corporation, 2002. 
NASA/TM—2014-218347 11 

5. Williams, K.V.; Vaziri, R.; and Poursartip, A.: “A physically based continuum damage mechanics 

model for thin laminated composite structures,” International Journal of Solids and Structures, Vol. 

40, pp. 2267-2300, 2003. 

6. Pinho, S.T.; Iannucci, L.; and Robinson, P.: “Physically-based failure models and criteria for 

laminated fibre-reinforced composites with emphasis on fibre kinking: Part I: Development,” 

Composites: Part A, Vol. 37, pp. 63–73, 2006. 

7. Pinho, S.T.; Iannucci, L.; and Robinson, P.: “Physically-based failure models and criteria for 

laminated fibre-reinforced composites with emphasis on fibre kinking: Part II: FE implementation,” 

Composites: Part A, Vol. 37, pp. 766–777, 2006. 

8. Maimi, P.; Camanho, P.P.; Mayugo, J.A.; and Davila, C.G.: “A continuum damage model for 

composite laminates: Part I – Constitutive model,” Mechanics of Materials, Vol. 39, pp. 897–908, 

2007. 

9. Maimi, P.; Camanho, P.P.; Mayugo, J.A.; and Davila, C.G.: “A continuum damage model for 

composite laminates: Part II – Computational implementation and validation,” Mechanics of 

Materials, Vol. 39, pp. 909-919, 2007. 

10. Daniel, I.M; and Ishai, O.: Engineering Mechanics of Composite Materials. Oxford University Press, 

New York, 2006. 

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

Composites,” Journal of Composite Materials, Vol. 23, pp. 1009-1020, 1989. 

12. Hoffarth, C.; Harrington, J.; Rajan, S.; Goldberg, R.; Carney, K.; DuBois, P.; Blankenhorn, G.: 

“Verification and Validation of a Three-Dimensional Generalized Composite Material Model,” 13th

International LS-DYNA Users Conference, Livermore Software Technology Corporation, 2014. 

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

14. Bednarcyk, B.A.; and Arnold, S.M.: 2002. “MAC/GMC 4.0 User’s Manual - Keywords Manual,” 

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

D.C., 2002.