Canio Hoffarth, Joseph Harrington, and Subramaniam D. Rajan

Arizona State University, Tempe, Arizona

Robert K. Goldberg and Kelly S. Carney

Glenn Research Center, Cleveland, Ohio

Paul DuBois

George Mason University, Fairfax, Virginia

Gunther Blankenhorn

Livermore Software Technology Corporation, Livermore, California

Veri¿ cation and Validation of a Three-Dimensional 

Generalized Composite Material Model

NASA/TM—2015-218446

January 2015
NASA STI Program . . . in Pro¿ le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scienti¿ 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 Of¿ 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 signi¿ cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signi¿ cant 

scienti¿ 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. Scienti¿ c 

and technical ¿ 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. Scienti¿ c and 

technical ¿ ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scienti¿ c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scienti¿ 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 scienti¿ c and 

technical material pertinent to NASA’s mission.

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services.

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

• Phone the NASA STI Information Desk at

 757-864-9658

• Write to:

NASA STI Information Desk

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199
Canio Hoffarth, Joseph Harrington, and Subramaniam D. Rajan

Arizona State University, Tempe, Arizona

Robert K. Goldberg and Kelly S. Carney

Glenn Research Center, Cleveland, Ohio

Paul DuBois

George Mason University, Fairfax, Virginia

Gunther Blankenhorn

Livermore Software Technology Corporation, Livermore, California

Veri¿ cation and Validation of a Three-Dimensional 

Generalized Composite Material Model

NASA/TM—2015-218446

January 2015

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

13th US and International LS-DYNA Users Conference

sponsored by the Livermore Software Technology Corporation

Dearborn, Michigan, June 8–10, 2014
Available from

NASA STI Information Desk

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

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

Trade names and trademarks are used in this report for identi¿ cation 

only. Their usage does not constitute an of¿ cial endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

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

Acknowledgments

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

Grant #12-G-001 entitled “Composite Material Model for Impact Analysis”, William Emmerling, Technical Monitor.
NASA/TM—2015-218446 1 

Verification and Validation of a Three-Dimensional Generalized 

Composite Material Model 

Canio Hoffarth, Joseph Harrington, and Subramaniam D. Rajan 

Arizona State University 

Tempe, Arizona 85004 

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 

Gunther Blankenhorn 

Livermore Software Technology Corporation 

Livermore, California 94551 

Abstract 

A general purpose orthotropic elasto-plastic computational constitutive material model has been 

developed to improve predictions of the response of composites subjected to high velocity impact. The 

three-dimensional orthotropic elasto-plastic composite material model is being implemented initially for 

solid elements in LS-DYNA as MAT213. In order to accurately represent the response of a composite, 

experimental stress-strain curves are utilized as input, allowing for a more general material model that can 

be used on a variety of composite applications. The theoretical details are discussed in a companion 

paper. This paper documents the implementation, verification and qualitative validation of the material 

model using the T800-F3900 fiber/resin composite material. 

Introduction 

Composite materials are now beginning to provide uses hitherto reserved for metals in structural 

systems such as airframes and engine containment systems, wraps for repair and rehabilitation and 

ballistic/blast mitigation systems. While material models exist that can be used to simulate the response of 

materials in these demanding structural systems, they are often designed for use with specific 

homogeneous materials such as metals (Refs. 1 and 2), polymers (Ref. 3) and wood (Ref. 4). While 

material models exist to simulate the impact response of composites, most are tailored specifically for a 

particular application and have limitations - purely elastic, no rate sensitivity, implementation for solid 

elements only, and limited damage and failure characterization (Refs. 5 to 10). 

A few examples of different approaches to composite damage modeling include Littell et al. 

(Ref. 11), Matzenmiller et al. (Ref. 12), and others (Refs. 13 to 18). Yen’s comprehensive model (Ref. 19) 

is implemented in LS-DYNA as *MAT_162. Availability of constitutive models for composite materials 

in LS-DYNA is discussed in some length in our companion paper (Ref. 20) and that explains the 

motivation for this study. 

 
NASA/TM—2015-218446 2 

In the next section we briefly discuss the theoretical development of the constitutive model. This is 

followed by details of the numerical algorithm, verification tests for a composite used in the aerospace 

industry, the use of the composite in a high speed contact example and finally, concluding remarks 

including details of our ongoing and future work. 

Orthotropic 3D Elasto-Plastic Composite Material Model 

The material law that the model is built upon describes the elastic and permanent deformation of the 

composite with full three-dimensional implementation suitable for solid elements. Only the relevant 

equations used in the numerical algorithm are presented here. The theoretical basis and discussions are 

available in a companion paper (Ref. 20). Current development of the model includes the complete rate 

independent elasto-plastic deformation model, with damage and failure to be added later. 

The commonly used Tsai-Wu composite failure model is generalized as a yield surface for the 

plasticity model and is defined as 

 

11 11 11 11 12 13

22 22 22 12 22 23

33 33 33 13 23 33

123

12 12 44 12

23 23 55 23

31 31 66 31

000

000

000

() 0 0 0 000 0 0

000 0 0

000 0 0

T FFF

FFF

FFF

f a FFF F

F

F

a o a o a oa o VV V

« » « » « »« » VV V

VV V V   « » « » « »«

VV V

VV V

¬ 1⁄4 ¬ 1⁄4 ¬ 1⁄4¬ 1⁄4 VV V

»

»

»

»

 (1) 

where a = –1 and 

1

11 11

2

22 22

3

33 33

1 1

1 1

1 1

T C

T C

T C

F

F

F

 

V V

 

V V

 

V V

 11

11 11

22

22 22

33

33 33

1

1

1

T C

T C

T C

F

F

F

 

V V

 

V V

 

V V

 44 2

12

55 2

23

66 2

31

1

1

1

F

F

F

 

V

 V

 

V

 (1a) 

 

2

45 45

2 1 ( ) , 1,2,3, 3

2

i j

ij ii jj kk i j i j

F F

F F F F ij k i  



     

V V

The stress components of the yield function coefficients correspond to the current yield stresses 

associated with a set of experimental stress-strain curves, which vary with the effective plastic strain, thus 

allowing the model to describe different hardening properties in each direction. The non-associative flow 

surface is defined as 

222 222

11 11 22 22 33 33 12 11 22 23 22 33 31 33 11 44 12 55 23 66 31 hH H H H H H H H H V  V  V  VV  V V  VV  V  V  V 222 (2) 

where the Hij terms are the constant flow rule coefficients. The rest of the details can be found in the 

companion paper (Ref. 20). 

The first six flow rule coefficients are computed directly from the assumed flow rule coefficient 

value, H11 and the plastic Poisson’s ratios. 
NASA/TM—2015-218446 3 

12 11 12

13 11 13

12

22 11

21

p

p

p

p

H H

H H

H H

 Q

 Q

Q 

Q

23 12

23 11

21

13

33 11

31

p p

p

p

p

H H

H H

Q Q Q

Q 

Q

 (3a) 

The last three flow rule coefficients (H44, H55, H66) are calculated iteratively with the user-defined 

range to best fit the input experimental shear stress versus shear strain curves (Ref. 20). For example, to 

determine H44, the following optimization problem is posed using the input curve for the “1-2” shear case. 

Find H44 to minimize 

  

2

44 22 12

1

( ) ˆ ˆ

n

i i

i

f H

 

a o

 V V « » ¬ 1⁄4 ¦ (3b) 

such that 

min max H HH 44 44 44 d d (3c) 

where n is the number of data points on the master curve,  22  ˆ i

V is the ith effective stress value from the 

master curve, and  12  ˆ i

V is the effective stress value for the shear curve using an assumed value of H44. 

Numerical Algorithm 

In this initial implementation, we do not account for rate and temperature dependence, damage 

accumulation nor failure. The following experimental data are needed as LS-DYNA input. 

(1) A total of 12 true stress versus true strain curves under quasi-static test conditions from (a) uniaxial 

tension in 1, 2 and 3-directions, (b) uniaxial compression in 1, 2 and 3-directions, (c) shear in 1-2, 

2-3 and 3-1 planes and (d) off-axis (e.g., 45q) uniaxial tension or compression in 1-2, 2-3 and 3-1 

planes. The tabulated x-y data for each curve is read in and used appropriately. 

(2) The modulus of elasticity, Poisson’s ratio and average plastic Poisson’s ratio obtained from the 

tension and compression tests. 

The effective plastic strain can be written in terms of the experimental stress versus total strain data. 

In the 1-direction for example, this is written as 

   

11

11 11 11

11 11

11 11 /

c c p

c pe

p p e

E

d h

§ · V 1⁄2

V H H   ̈  ̧°

° © 13⁄4 V H

°

H V H 3 °¿

 (4) 

where 11

c V is the experimental compressive true stress in the 1-direction, 11 H is the total true strain in the 

1-direction, E11 is the elastic modulus in the 1-direction, 11

p H is the true plastic strain in the 1-direction, 
NASA/TM—2015-218446 4 

p

e H is the effective plastic strain and h is the current value of the plastic potential function, which is 

equivalent to the effective stress, shown in Equation (2). The material model uses a typical elastic stress 

update as 

  1 : nn p t  ı ı '  C İ İ   (5) 

where C is the orthotropic stiffness matrix, 't is the time step, İ is the total strain rate and p İ is the 

plastic strain rate. The stiffness matrix is written in terms of the compliance matrix as 

1

21 31

11 22 33

32

22 33

1 33

12

23

31

1 000

1 000

1 000

1 0 0

1 0

1

v v

EEE

v

E E

E

G

Sym G

G





a o   « »



 

¬ 1⁄4

C S (6) 

where Eij are the elastic moduli in the principal material directions, Gij are the elastic shear moduli and 

ij v are the elastic Poisson’s ratios. The yield stresses used to determine the flow rule coefficients are 

summarized into a single vector, corresponding to each of the experimental test curves as 

11 22 33 11 22 33 12 23 31 45 12 45 23 45 31

T TT T CC C C C C  V V V V V V V V V V V V a o ¬ 1⁄4 q . Lastly, the consistency 

condition is written in terms of the rate of yield function as follows 

0 f f f w w  

w w

ı q

ı q

   (7) 

where the terms are as defined earlier. Equation (7) can be expanded as follows by introducing the 

elasto-plastic constitutive law 

:: 0 f h fd f d

w ww § · O O  ̈  ̧ w wwO © 1

q C İ C

ı ı q

   (8) 

where ı is written in terms of the constitutive matrix and strain rates and O is the effective plastic strain. 

Solving for the effective plastic strain rate produces the following equation. 
NASA/TM—2015-218446 5 

:

:

f

f h fd

d

w

w O w ww



w w wO

C İ

ı

q C

ı ı q

  (9) 

To numerically solve for the increment in effective plastic strain, a radial return algorithm is used. 

Initially, a perfectly elastic response is assumed. Thus, an elastic predictor is used to compute the stresses 

as 

        

        

       

   

   

   

11 11 11 11 12 22 13 33

22 22 12 11 22 22 23 33

33 33 13 11 23 22 33 33

12 12 44 12

23 23 55 23

31 31 66 31

e n

e n

e n

e n

e n

e n

Ct Ct Ct

Ct Ct Ct

Ct Ct Ct

C t

C t

C t

V V  ' H  ' H  ' H

V V  ' H  ' H  ' H

V V  ' H  ' H  ' H

V V  ' H

V V  ' H

V V  ' H

 



 







 (10) 

With the elastic trial stresses computed, a trial yield function value can be calculated and used to 

determine if the load step is elastic or plastic. It should be noted that the elastic trial stress is used as the 

stress in the next step if the trial yield function value is less than or equal to zero. Otherwise, the stress 

state is beyond the yield surface and must be brought back using a radial return and plastic corrector, 'O . 

If the trial yield function happens to be greater than zero, then 'O must be greater than zero. The value of 

'O is determined using a secant iteration, with 1 'O 0 for the first iteration and the second value for the 

increment of the plastic multiplier determined from consistency, 

    2

: :

: :

en en

e e

ee ee

f f

f h fd f h

d

w w   w w 'O | w ww w w



w w wO w w

ıı ıı

ı ı

q C C

ı ı q ı ı

 (11) 

where the derivatives of q are taken as zero, meaning that the response is perfectly plastic. The initial 

calculation of the second increment corresponds to perfect plasticity in order to ensure that the stress state 

returns to the interior of the yield surface, resulting in a negative value of the yield function and bounding 

the solution. From each estimate of 'O , the corresponding corrected stresses can be computed as well as 

the yield stresses obtained from the input stress versus effective plastic strain (O) data. With these terms 

calculated for both estimates of 'O , the yield function value for each can be computed with the third 

estimate evaluated as 

2 1

3 11

2 1 f

f f

'O  'O 'O 'O   (12) 

The corrected stresses and yield stresses associated with ǻȜ3

 are then calculated and used to obtain a 

new estimate of the yield function value. Convergence of the secant iteration is determined by the 

following conditions, where f3 is the value of the yield function given the effective plastic strain ǻȜ3

, 
NASA/TM—2015-218446 6 

13 11

3

333 22 23 fff 0 ;0 ;0

­ ­ ° ° 'O 'O 'O 'O

!    |  'O 'O ® ®

° °  ̄  ̄ 'O 'O 'O 'O

 (13) 

with the secant iteration continuing if the new yield function value is not less than some defined tolerance. 

Once convergence is met, the plastic multiplier increment to return the stress state to the yield surface is 

known and the stresses can be updated as 

1 : n e

e

 wh  'O

w

ı ı C

ı

 (14) 

Finally, the yield stresses are updated as well, using the new value of the overall effective plastic 

strain, O, in each input curve to determine the corresponding yield stress level as 

  n n 1 3 q q  O  'O (15) 

This results in anisotropic strain hardening, as a yield stress increase in each direction is initiated with 

plasticity in any direction, but at different levels. A detailed algorithm is presented below. 

Step 1: Preprocessing

(a) Convert stress/strain input curves to stress versus effective plastic strain using Equation (4). 

(b) Store initial yield stresses in q . 

(c) Calculate and store the elastic moduli from the initial yield stress and strain values. 

(d) Store the three elastic Poisson’s ratio values and flow rule coefficients from input. 

(e) Compute optimal values of the flow rule coefficients so as to match the input curves as closely as 

possible. 

Step 2: Initialization

Parameters available at start of step: qn , ın , (, ) nn n İ İ  ' o' t , On . 

Step 3: Elastic predictor 

Compute and set the yield function coefficients using Equation (1a). 

Construct the constitutive matrix using Equation (6). 

Compute elastic trial stresses, 1 e ın , using Equation (10). 

Compute the trial yield function, 1

trial

nf  , using the elastic trial stresses in Equation (1). 

If 1 0 trial

nf  d , the elastic solution is correct, set 0 'O n and go to stress update. Else go to plastic 

corrector. 

Step 4: Plastic corrector

Set 1 'O 0 . 

Calculate 2 'O from Equation (11). 

Loop through secant iteration for specified max number of iterations: 

Compute the new estimate of the stress for each plastic multiplier increment   1 2 'O 'O , using 

Equation (14). 
NASA/TM—2015-218446 7 

Calculate the effective plastic strains   1 2 O O, at the next time step using the estimates of the plastic 

multiplier increment; then obtain the strain hardening parameters with the effective plastic strain 

from the curve data (true stress as a function of effective plastic strain). 

Compute new plastic multiplier increment, 3 'O . Calculate new estimate for the yield function, 3 f . 

If 3 f d tol , set 3 'O 'O , exit the loop and go to stress update. Else update secant iteration 

parameters using Equation (14) and go to next step of secant iteration. 

Step 5: Stress Update

Calculate ın1 using Equation (14). 

Step 6: Update history variables for plastic work and work hardening parameters ( q and O ) 

Set O O  'O n nn 1 . 

Determine new yield stresses, qn1, using On1 and Equation (15). 

Model Verification 

The composite material model was tested and verified using experimental data obtained from 

T800S/3900-2B [P2352W-19] BMS8-276 Rev-H-Unitape fiber/resin unidirectional composite (Ref. 21). 

Toray describes T800S as an intermediate modulus, high tensile strength graphite fiber. The epoxy resin 

system is labeled F3900 where a toughened epoxy is combined with small elastomeric particles to form a 

compliant interface or interleaf between fiber plies to resist impact damage and delamination (Ref. 22). 

Magnified views of the composite are shown in Figures 1 and 2. 

Figure 1.—Side view obtained 

using optical microscopy. 

Figure 2.—Logitudinal view obtained 

using a SEM. 

 
NASA/TM—2015-218446 8 

The source of the input data for the model (both experimental and virtual) is listed in Table 1. In the 

Data column, Experimental refers to experimental data generated at Wichita State University (Ref. 21). 

For cases where actual experimental data were not available, numerical simulations (virtual experiments) 

were conducted to obtain the required stress-strain curves. To conduct the numerical experiments, 

micromechanics analyses were conducted where the properties of the fiber and matrix were used to 

simulate the response of the overall composite. Two methods were used to conduct the numerical 

experiments. First, the NASA GRC developed code MAC-GMC was used (Ref. 23). MAC-GMC utilized 

the Generalized Method of Cells to conduct a fully analytical simulation of the composite response. A 

second approach used a Virtual Testing Software System (Ref. 24) in which finite element models were 

constructed and analyzed where the fiber and matrix were explicitly modeled. For this particular 

unidirectional composite, as indicated in the table, due to the transverse isotropy of the material the 

3-direction response was set equal to the 2-direction response. Likewise, the 1-3 shear response was set 

equal to the 1-2 shear response, and the off-axis response in the 1-3 plane was set equal to the off-axis 

response in the 1-2 plane. More details on the procedures used in the numerical experiments will be given 

in a future report. 

The fiber (transversely isotropic, linear elastic) and the matrix (isotropic, elastic perfectly plastic) 

properties used in the numerical experiments are listed in Table 2. The properties for the latter are not 

publicly available. Hence they were correlated to obtain the values shown in the table. 

The values of the coefficients in the flow rule were determined from the experimental (actual or 

virtual) data and the procedures outlined in Equations (3a), (3b), and (3c). The correlation curves that 

were used to optimize the shear coefficients in the flow law are shown in Figures 3 and 4. In the plots, the 

[90°] transverse tensile curve was chosen as the baseline master curve for the optimization. 

TABLE 1.—GENERATION OF INPUT DATA FOR 

T800-F3900 COMPOSITE 

Curve Data 

Tension Test (1-Direction) Experimental 

Tension Test (2-Direction) MAC-GMC and VTSS 

Tension Test (3-Direction) Transverse isotropy 

Compression Test (1-Direction) Experimental 

Compression Test (2-Direction) MAC-GMC and VTSS 

Compression Test (3-Direction) Transverse isotropy 

Pure Shear Test (1-2 Plane) Experimental 

Pure Shear Test (2-3 Plane) MAC-GMC and VTSS 

Pure Shear Test (1-3 Plane) Transverse isotropy 

Off-Axis Test (45°, 1-2 Plane) MAC-GMC and VTSS 

Off-Axis Test (45°, 2-3 Plane) MAC-GMC and VTSS 

Off-Axis Test (45°, 1-3 Plane) Transverse isotropy 

TABLE 2.—T800-F3900 FIBER AND MATRIX PROPERTIES 

(VOLUME FRACTION = 0.54) 

Property Fiber Matrix 

E1, psi 4(107) 5(105) 

E2, psi 2.25(107) 

E3, psi 2.25(107) 

Q12 0.2 0.35 

Q23 0.25 

Q13 0.25 

G, psi 1.5(107) 1.85(105) 

Vy, psi ---------- 2(104) 
NASA/TM—2015-218446 9 

Figure 3.—Comparison of experimental curve with optimized H44 (and H66) value curve. 

Figure 4.—Comparison of experimental curve with optimized H55 value curve. 

 
NASA/TM—2015-218446 10 

After the determination of the required constants, the verification process was executed. In the 

verification test results shown in the figures that follow, the term Experimental (Input) refers to 

experimental data (when available) or to data created via virtual testing, and the term Simulated 

(x-element) refers to the data created using results from finite element simulations where x is the number 

of finite elements in the FE model. 

Schematics for the tension and compression tests are shown in Figures 5 and 6, respectively. The 

simulations for the tension and compression tests in all three directions were executed using different 

mesh sizes ranging from 1 to 64 elements. The 64-element meshes for the tension and compression test 

cases are shown in Figures 7 and 8, respectively. 

2”

8”

A t=0.25”

B

C

D

E

1

2

 (a) 

2”

8”

A t=0.25”

B

C

D

E

2, 3

1

(b) 

Figure 5.—Schematics for tension test simulation cases (a) 1-direction (b) 2 and 3-directions. 

1 (x)

0.25”

2 (y)

3 (z)

0.25”

0.25”

(a) 

2/3 (x)

0.25”

1 (y)

-3/2 (z)

0.25”

0.25”

(b) 

Figure 6.—Schematics for compression test simulation cases (a) 1-direction (b) 2 and 3-directions. 

Figure 7.—The 64-element mesh for tension test cases. Figure 8.—The 64-element mesh 

for compression test cases. 

 
NASA/TM—2015-218446 11 

For both tension and compression tests, nodes on face ABC were constrained in the x-direction and 

the center node was also constrained in the y-direction. A 0.5 in./s displacement in the x-direction was 

applied to all the nodes on the right face, E-D. For the compression tests, care was taken to ensure that the 

model yielded but did not buckle. The simulated and experimental stress-strain curves for the tension and 

compression tests are shown in Figures 9 and 10. In all cases the correlation between the input curves and 

the model predictions was excellent. Note that the transverse tension and compression curves (both the 

numerical experiments and model simulations) were artificially extended well into the nonlinear range. 

The extension was required in order to both fully characterize the model and to explore the capability of 

the model to correctly simulate the nonlinear response of the composite. 

Figure 9.—Simulated and experimental stress-strain curves for 1-direction (a) tension and (b) compression. 
NASA/TM—2015-218446 12 

Figure 10.—Simulated and experimental stress-strain curves for 2/3-directions (a) tension and (b) compression. 

 
NASA/TM—2015-218446 13 

A schematic for the pure shear test is shown in Figure 11. Nodes on the bottom surface are restrained 

in the x and y-directions while nodes on the top surface are restrained in the y-direction. A displacement 

controlled loading of 0.5 in./s in the x-direction is applied to all the nodes on the top surface. The 

64-element mesh is shown in Figure 12. The simulated and experimental stress-strain curves in the 

1-2/3-1 plane are shown in Figure 13 and in the 2-3 plane in Figure 14. Once again, the correlation 

between the input curves and the model predictions was excellent. 

The (tension) off-axis test model in the 1-2 plane is shown in Figure 15. Nodes on face ABC were 

constrained in the x-direction and the center node was also constrained in the y-direction. The local 

element coordinate system was rotated 45° in order to simulate the off-axis response. The simulated and 

experimental stress-strain curves are compared in Figure 16, showing very little differences in the results, 

which at least partially demonstrates the ability of the model simulate the response of a composite when it 

is subjected to a multi-axial stress state (in the material axis system). Finally, Figure 17 shows the results 

from a loading-unloading-reloading test case. The model captured the results expected for a pure 

plasticity formulation, including a linear elastic unload with the original modulus, and a resumption of 

plastic flow upon reloading to the original unload point. Capturing the reduction in the unloading 

modulus expected for materials of this type will require the addition of a complementary damage model 

to the formulation. 

1 (x)

2”

2”

0.25”

2 (y)

3 (z)

Figure 11.—Schematic for pure shear test in 

1-2/3-1 plane. 

Figure 12.—The 64-element mesh for pure shear 

test cases.

Figure 13.—Simulated and experimental stress-strain curves for pure shear in the 1-2/3-1 plane. 
NASA/TM—2015-218446 14 

Figure 14.—Simulated and experimental stress-strain curves for pure shear in the 2-3 plane. 

2”

8”

A t=0.25”

B

C

D

E

1

2

x

y

Figure 15.—Schematic for 45° off axis test in 1-2/3-1 plane. 

Figure 16.—Simulated and experimental stress-strain curves for 45° off axis test in 1-2/3-1 plane. 
NASA/TM—2015-218446 15 

Figure 17.—Simulated and experimental stress-strain curves for unloading/reloading in the 2-direction. 

In summary, simulations of the curves represented by the 12 input cases as an initial verification 

study shows very good correlation with the input data. We have used boundary conditions that resemble 

the experimental test setup as closely as possible so as to capture the actual test conditions. The use of 

four different meshes with increasing mesh density, shows that the results are mesh independent. 

Structural Test 

To gage the ability of the model to produce reasonable results under more realistic structural loading 

conditions, a structural test case was designed to qualitatively validate the developed model. A simplysupported 12 in. long by 2.4 in. wide by 0.6 in. thick plate made of T800-F3900 composite material is 

subjected to an impact from a steel ball. Figure 18 shows the schematic diagram. The nodes on the bottom 

of the left edge are constrained in the x, y and z-directions while nodes on the bottom of the right edge 

were constrained only in the y-direction. In the absence of experimental data, the primary objectives were 

to ensure that the material model (in the absence of damage and failure) would yield realistic results. That 

is, the ball would impact the plate and bounce back, and the peak stresses in the plate would occur around 

the point of contact immediately after the impact. The LS-DYNA FE model is shown in Figure 19. 

References to LS-DYNA keywords can be found in Hallquist (Ref. 25). 

The model has 6,478 nodes and 5,707 elements of which 1,147 8-noded hexahedral elements are used 

to model the plate (master surface) and the remaining 4,560 elements are used in the ball (slave surface). 

The ball is modeled using *MAT_3 and the plate as *MAT_213. *CONTACT_AUTOMATIC_ 

SURFACE_TO_SURFACE is used to describe the contact between the ball and plate. Both the static (fs) 

and the dynamic (fd) coefficients of friction are set to 0.1. A 2 percent viscous damping coefficient is 

used. Segment-based contact is used by setting the soft parameter to 2. The steel ball is given a –700 in./s 

initial velocity in the y-direction. The analysis is carried out for 0.1 sec. The results are as expected. The 

exaggerated deformed shape is shown in Figure 20. The von Mises stress throughout the test duration is 

plotted for the elements closest to the point of contact in Figure 21. 
NASA/TM—2015-218446 16 

1.20"

SIDE VIEW

PLAN VIEW

12.00"

0.60" 2.40"

1.00"

v0

Figure 18.—Schematic diagram of the test. 

Figure 19.—LS-DYNA FE model at the start of the 

analysis. 

Figure 20.—Transient shape after impact. 

Figure 21.—von Mises stress versus time for elements near point of impact. 
NASA/TM—2015-218446 17 

Conclusions 

A generalized three-dimensional composite material model is developed and verified. The preliminary 

results are encouraging as the implemented constitutive model is able to reproduce the experimental stressstrain curves. In addition, a simple contact problem involving the composite material yields reasonable 

results. Future work includes the following: adding damage accumulation and failure in the material model, 

supporting temperature and strain rate dependent stress-strain curves, using alternate methods to compute 

the plastic multiplier when the initial yield and flow surface convexity conditions are not met, using a fiberreferenced convected system to track the principal material directions and supporting plate/shell elements 

and to validate the developed model using data from several structural tests. 

References 

1. L. P. Moreira and G. Ferron, “Finite element implementation of an orthotropic plasticity model 

for sheet metal forming simulations,” Latin American Journal of Solids and Structures, vol. 4, 

pp. 149-176, 2007. 

2. M. Ganjiani, R. Naghdabadi and M. Asghari, “An elastoplastic damage-induced anisotropic 

constitutive model at finite strains,” International Journal of Damage Mechanics, vol. 22, 

pp. 499-529, 2012. 

3. J. S. Bergstrom, “Constitutive Modeling of Elastomers - Accuracy of Predictions and Numerical 

Efficiency,” 2005. [Online]. Available: PolymerFEM.com. 

4. A. Tabiei and J. Wu, “Three-dimensional nonlinear orthotropic finite element material model for 

wood,” Composite Structures, vol. 50, pp. 143-149, 2000. 

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

Composites,” Journal of Composite Materials, 1989. 

6. R. Vaziri, M. D. Olson and D. L. Anderson, “A Plasticity-Based Constitutive Model for 

Fibre-Reinforced Composite Laminates,” Journal of Composite Materials, pp. 512-535, 1991. 

7. G. A. Holzapfel and T. C. Gasser, “A viscoelastic model for fiber-reinforced composites at finite 

strains: Continuum basis, computational aspects and applications,” Computer methods in applied 

mechanics and engineering, vol. 190, pp. 4379-4403, 2001. 

8. P. B. Lourenco, R. D. Borst and J. G. Rots, “A Plane Stress Softening Plasticity Model for 

Orthotropic Materials,” International Journal for Numerical Methods in Engineering, pp. 4033-4057, 

1997. 

9. Y. D. Murray and L. E. Schwer, “Verification of a general-purpose laminated composite shell 

element implementation: Comparisons with analytical and experimental results,” Finite Elements in 

Analysis and Design, vol. 13, no. 1, pp. 1-16, 1993. 

10. G. C. Saha and A. L. Kalamkarov, “Micromechanical Thermoelastic Model for Sandwhich 

Composite Shells made of Generally Orthotropic Materials,” Journal of Sandwich Structures and 

Materials, vol. 11, pp. 27-56, 2009. 

11. J. D. Littell, W. K. Binienda, W. A. Arnold, G. D. Roberts and R. K. Goldberg, “Effect of 

Microscopic Damage Events on Static and Ballistic Impact Strength of Triaxial Braid Composites,” 

NASA, Washington, 2010. 

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

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

13. X. Xiao, “Modeling Energy Absorbtion with a Damage Mechanics Based Composite Material 

Model,” Journal of Composite Materials, vol. 43, pp. 427-444, 2009. 

14. F. Wu and W. Yao, “A fatigue damage model of composite materials,” International Journal of 

Fatigue, pp. 134-138, 2010. 

15. H.-S. Chen and S.-F. Hwang, “A Fatigue Damage Model for Composite Materials,” Polymer 

Composites, pp. 301-308, 2009. 
NASA/TM—2015-218446 18 

16. B. A. Gama, J.-R. Xiao, M. J. Haque, C.-F. Yen and J. W. J. Gillespie, “Experimental and 

Numerical Investigations on Damage and Delamination in Thick Plain Weave S-2 Glass Composites 

Under Quasi-Static Punch Shear Loading,” U.S. Army Research Laboratory, Aberdeen Proving 

Ground, 2004. 

17. A. F. Johnson, A. K. Pickett and P. Rozycki, “Computational methods for predicting impact damage 

in composite structures,” Composites Science and Technology, vol. 61, pp. 2183-2192, 2001. 

18. J. Ozbolt, V. Lackovic and J. Krolo, “Modeling fracture of fiber reinforced polymer,” International 

Journal of Fracture, vol. 170, pp. 13-26, 2011. 

19. C.-F. Yen, “Ballistic Impact Modeling of Composite Materials,” in 7th International LS-DYNA 

Users Conference, 2002. 

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

“Theoretical Development of an Orthotropic Three-Dimensional Elasto-Plastic Generalized 

Composite Material Model,” in 13th International LS-DYNA User’s Conference, Detroit, MI, 2014. 

21. K. S. Raju and J. F. Acosta, “Crashworthiness of Composite Fuselage Structures - Material Dynamic 

Properties, Phase I,” U.S. Department of Transportation: Federal Aviation Administration, 

Washington, DC, 2010. 

22. D. L. Smith and M. B. Dow, “Properties of Three Graphite/Toughened Resin Composites,” 

NASA Technical Paper 3102, September 1991. 

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

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

Washington, D.C., 2002. 

24. J. Harrington and S. D. Rajan, Test Results from Virtual Testing Software System, Technical 

Report, School of Sustainable Engineering and the Built Environment, ASU, March 2014. 

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

Corporation, Livermore, CA, 2013.