Autodyn Composite Modeling Guide

ANSYS, Inc. Release 2020 R1

Southpointe January 2020
2600 ANSYS Drive
Canonsburg, PA 15317 ANSYS, Inc. and
ansysinfo@[Link] ANSYS Europe,
Ltd. are UL
[Link] registered ISO
(T) 724-746-3304 9001: 2015
(F) 724-514-9494 companies.
Copyright and Trademark Information

© 2020 ANSYS, Inc. Unauthorized use, distribution or duplication is prohibited.

ANSYS, ANSYS Workbench, AUTODYN, CFX, FLUENT and any and all ANSYS, Inc. brand, product, service and feature
names, logos and slogans are registered trademarks or trademarks of ANSYS, Inc. or its subsidiaries located in the
United States or other countries. ICEM CFD is a trademark used by ANSYS, Inc. under license. CFX is a trademark
of Sony Corporation in Japan. All other brand, product, service and feature names or trademarks are the property
of their respective owners. FLEXlm and FLEXnet are trademarks of Flexera Software LLC.

Disclaimer Notice

THIS ANSYS SOFTWARE PRODUCT AND PROGRAM DOCUMENTATION INCLUDE TRADE SECRETS AND ARE CONFID-
ENTIAL AND PROPRIETARY PRODUCTS OF ANSYS, INC., ITS SUBSIDIARIES, OR LICENSORS. The software products
and documentation are furnished by ANSYS, Inc., its subsidiaries, or affiliates under a software license agreement
that contains provisions concerning non-disclosure, copying, length and nature of use, compliance with exporting
laws, warranties, disclaimers, limitations of liability, and remedies, and other provisions. The software products
and documentation may be used, disclosed, transferred, or copied only in accordance with the terms and conditions
of that software license agreement.

ANSYS, Inc. and ANSYS Europe, Ltd. are UL registered ISO 9001: 2015 companies.

U.S. Government Rights

For U.S. Government users, except as specifically granted by the ANSYS, Inc. software license agreement, the use,
duplication, or disclosure by the United States Government is subject to restrictions stated in the ANSYS, Inc.
software license agreement and FAR 12.212 (for non-DOD licenses).

Third-Party Software

See the legal information in the product help files for the complete Legal Notice for ANSYS proprietary software
and third-party software. If you are unable to access the Legal Notice, contact ANSYS, Inc.

Published in the U.S.A.

Table of Contents
1. Composite Material Modeling Introduction ........................................................................................... 1
2. Principal Directions in Autodyn .............................................................................................................. 3
2.1. Autodyn-2D ...................................................................................................................................... 3
2.1.1. X-Y Space ................................................................................................................................. 3
2.1.2. Polar Space .............................................................................................................................. 4
2.1.3. I-J-K Space ............................................................................................................................... 4
2.2. Autodyn-3D Lagrange and ALE Parts ................................................................................................. 5
2.2.1. X-Y-Z Space ............................................................................................................................. 5
2.2.2. I-J-K Space ............................................................................................................................... 6
2.3. Autodyn-3D Shell Parts ..................................................................................................................... 7
2.3.1. X-Y-Z Space ............................................................................................................................. 9
2.3.2. I-J-K Space ............................................................................................................................... 9
3. Orthotropic Constitutive Models .......................................................................................................... 11
3.1. Orthotropic Elastic Model ............................................................................................................... 11
3.2. Equations of State ........................................................................................................................... 14
3.2.1. Shock Equation of State .......................................................................................................... 16
3.2.2. Polynomial Equation of State .................................................................................................. 16
3.3. Strength Models ............................................................................................................................. 17
3.3.1. Elastic .................................................................................................................................... 17
3.3.2. VonMises ............................................................................................................................... 17
3.3.3. Orthotropic Yield Strength and Hardening Model ................................................................... 17
[Link]. Quadratic Limit Surface ................................................................................................. 18
[Link]. Stress Return Algorithm ................................................................................................. 20
4. Orthotropic Material Failure Models .................................................................................................... 23
4.1. Brittle Damage Model ..................................................................................................................... 23
4.1.1. Failure Initiation ..................................................................................................................... 23
4.1.2. Post-Failure Response ............................................................................................................ 24
[Link]. Isotropic Post-Failure ..................................................................................................... 24
[Link]. Orthotropic Post-Failure Response ................................................................................. 25
4.2. Orthotropic Damage Model ............................................................................................................ 28
4.2.1. Failure Initiation ..................................................................................................................... 28
4.2.2. Damage Model ...................................................................................................................... 28
5. Material Characterization Tests ............................................................................................................ 31
5.1. Directional Strength Properties ....................................................................................................... 31
5.1.1. In-Plane Tension Tests ............................................................................................................. 31
5.1.2. Out-of-Plane Shear ................................................................................................................. 33
5.2. Equation of State Properties: Inverse Flyer Plate Tests ....................................................................... 34
5.3. Delamination Properties ................................................................................................................. 35
5.3.1. Direct Plate Impact - Mode I Delamination Strength ................................................................ 36
5.3.2. Double Cantilever Beam – Mode I Delamination Energy .......................................................... 37
5.3.3. Double Notch Shear – Mode II Delamination Strength ............................................................. 38
5.3.4. End Notched Flexure - Mode II Delamination Energy ............................................................... 40
6. Derivation of Material Properties ......................................................................................................... 43
6.1. Constitutive Properties ................................................................................................................... 43
6.1.1. Elastic Properties .................................................................................................................... 43
6.1.2. Plasticity Parameters .............................................................................................................. 45
6.2. Equation of State Parameters .......................................................................................................... 47
6.3. Failure and Softening Properties ..................................................................................................... 49
6.3.1. Failure Properties ................................................................................................................... 49
6.3.2. Softening Properties .............................................................................................................. 49

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. iii
Composite Modeling Guide

6.4. Calculating Laminate Properties from Uni-Directional Data .............................................................. 50

7. Example Applications ........................................................................................................................... 53
7.1. Hypervelocity Impacts .................................................................................................................... 53
7.1.1. Advanced Material Model for Hypervelocity Impact Simulations .............................................. 53
7.1.2. Advanced Material Damage Models ........................................................................................ 55
7.2. Ballistic Impact Examples ................................................................................................................ 58
7.2.1. Ballistic Fragment Impacts on Aramid Composite Plates .......................................................... 59
7.2.2. Polyethylene Fiber-Based Armor ............................................................................................. 59
7.3. Bird Strike Example using Composite Shell Elements ....................................................................... 60
8. Recommendations ................................................................................................................................ 63
9. References ............................................................................................................................................. 65

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
iv of ANSYS, Inc. and its subsidiaries and affiliates.
List of Figures
2.1. Material directions in X-Y space ............................................................................................................... 3
2.2. Material directions in polar space ............................................................................................................ 4
2.3. Material directions in I-J space ................................................................................................................. 4
2.4. Material directions in Autodyn-3D defined in XYZ-space .......................................................................... 5
2.5. Material directions in Autodyn-3D defined in XYZ-space .......................................................................... 6
2.6. Material directions in Autodyn-3D defined in IJK-space ............................................................................ 7
2.7. Co-rotational coordinate system used for 3D structured Shell elements ................................................... 8
2.8. Co-rotational coordinate system used for 3D unstructured Shell elements ............................................... 8
2.9. Material directions defined in XYZ-space for 3D Shell elements ................................................................ 9
2.10. Material directions defined in IJK-space for 3D structured Shell .............................................................. 9
3.1. Typical In-Plane Stress-Strain Behavior of Kevlar-Epoxy .......................................................................... 11
3.2. Computational Cycle ............................................................................................................................. 20
3.3. Yield surface ......................................................................................................................................... 21
3.4. Schematic representation of backward-Euler return algorithm ............................................................... 22
4.1. Schematic Illustration of Crack Softening Algorithm .............................................................................. 29
5.1. Schematic Diagram of a Tensile Test Specimen ...................................................................................... 32
5.2. Calculation of In Plane and Out of Plane Poisson Ratios .......................................................................... 32
5.3. Typical Recorded Uniaxial Stress-Strain Relationship .............................................................................. 33
5.4. Short Beam Shear Test Setup ................................................................................................................. 34
5.5. Typical Shear Stress - Shear Angle Relationship ...................................................................................... 34
5.6. Configuration of Inverse Planar Impact Experiments and a Typical Velocity Trace from the Rear Surface
of the Witness Plate .................................................................................................................................... 35
5.7. Delamination Modes. (a) Mode I : Normal Delamination. (b) Mode II : Shear Delamination. (c) Mode III :
Shear Delamination .................................................................................................................................... 36
5.8. Experimental Configuration of the Direct Plate Impact Experiment ........................................................ 36
5.9. Typical Velocity Trace from a Direct Plate Impact Test ............................................................................. 37
5.10. Illustration of the Double Cantilever Beam Test .................................................................................... 38
5.11. Force – Displacement Curve for Determination of the Mode I Fracture Energy Release Rate .................. 38
5.12. Configuration of short beam shear test sample .................................................................................... 39
5.13. Typical Output from DNS Test: Interlaminar Shear Strength .................................................................. 39
5.14. Schematic of the ENF Configuration Used to Determine the Interlaminar Fracture Energy GIIC .............. 40
5.15. Typical Load Displacement Curve for Determination of the Mode II Interlaminar Fracture Energy .......... 41
6.1. Derivation of Master Relationship from Uniaxial Tension Test ....................................................... 46
6.2. Longitudinal Versus Transverse Strain Measured in 0° Tension Tests. Red Line Indicates Region used to
Calculate Value for In-Plane Poissons Ratio. (Picture Courtesy of EMI [2]) ....................................................... 46
6.3. Results from Simulations of 45° Tension Tests. (Experimental Result Courtesy of EMI [2]) ......................... 47
6.4. Kevlar/Epoxy IFPT, Influence of Shock Effects. (Experimental Result Courtesy of EMI [1]) .......................... 48
6.5. Shock Velocity versus Particle Velocity Relationship for Kevlar-epoxy Inverse Flyer Plate Tests. (Data
courtesy of EMI [1]) ..................................................................................................................................... 48
7.1. Inverse Flyer Plate Tests, Experimental and AMMHIS Model Results ......................................................... 53
7.2. Alenia/EMI Test A8611 – Material Status During Impact on Reference Shielding, 15mm Diameter Projectile,
6.5km/s ...................................................................................................................................................... 54
7.3. Alenia/EMI Test A8611- Key Features of Material Response During Impact on Reference Shielding, 15mm
Diameter Projectile, 6.5km/s ........................................................................................................................ 54
7.4. Results from Simulation of 0° Tension Test ............................................................................................. 55
7.5. Simulation of Short Beam Shear Test ..................................................................................................... 55
7.6. Schematic Diagram of Experimental Setup ............................................................................................ 56
7.7. Simulation Results of 276 m/s Impact Velocity ....................................................................................... 56
7.8. Test 4355 – Through Thickness Damage During Impact .......................................................................... 57
7.9. Final Damage of Test 4355 ..................................................................................................................... 58

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. v
Composite Modeling Guide

7.10. Fragment Impact on KFRP at 483m/s, Simulation and Experimental Results .......................................... 59
7.11. Autodyn Simulations of 1.1g FSP Impacting 3.2mm Dyneema UDHB25 (Magenta Regions Indicate
Delamination) ............................................................................................................................................ 60
7.12. Initial Configuration of Composite Tail Section ..................................................................................... 61
7.13. Initial Configuration of Composite Tail Section: Impact Zone ................................................................ 61
7.14. Material Status Following Bird Strike .................................................................................................... 61

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
vi of ANSYS, Inc. and its subsidiaries and affiliates.
List of Tables
4.1. Orthotropic Post-Failure Options ........................................................................................................... 25
6.1. Derivation of Orthotropic Material Elastic Properties .............................................................................. 44
6.2. Failure Properties and the Tests Used to Measured Them ....................................................................... 49
6.3. Fracture Energies and Tests From Which They May be Derived ............................................................... 49

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. vii
 Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
viii of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 1: Composite Material Modeling Introduction
The Autodyn hydrocode has extensive capabilities for the modeling of composite materials subjected
to a range of loading conditions. A simple linear-elastic orthotropic constitutive model, inherent in
which is a linear equation of state, suitable for modeling applications subjected to structural (rather
than shock) type loading can be used. Or, for applications such as hypervelocity impacts where the
shock effects are obviously important, the orthotropic model can be coupled with nonlinear equations
of state. Damage/Failure can be treated as brittle via directional failure models like Material Stress and/or
Strain Failure (see Brittle Damage Model (p. 23)). It can also be treated by a specific Orthotropic Damage
model (see Orthotropic Damage Model (p. 28)), which takes into account the softening behavior that
can sometimes be observed when failure occurs in composites.

This document provides a single point of reference for using these models and describes all of the
above options. Often the greatest difficulty in using such models is a lack of available material data or
not knowing how to properly characterize the material should experimental facilities be available.
Therefore this document also provides a description of material characterization experiments that may
be performed in order to calculate the required material input parameters.

ANSYS Autodyn is an integrated analysis program designed for nonlinear dynamics problems. There is
at times a need to consider the treatment of materials where the properties of materials are not
identical in all directions (for example, composite laminates, fiber reinforced materials). Material models
suitable for such anisotropic material behavior have been developed in Autodyn. These models have
differing levels of complexity and require differing amounts of material data as input.

The purpose of this document is to fully describe the composite material modelling capabilities of
Autodyn. Further guidance is offered to help the user select the most appropriate models for a specific
application and to obtain the relevant material data.

After describing the definition of the principal directions in Autodyn in Principal Directions in Auto-
dyn (p. 3), the material models are presented for both the constitutive behavior and for predicting
failure/damage in Orthotropic Constitutive Models (p. 11) and Orthotropic Material Failure Models (p. 23)
respectively. The most common problems in using such models are the lack of material data and the
difficulty in obtaining the measurements. In Material Characterization Tests (p. 31) a series of experi-
mental tests are presented which have been used in projects at ANSYS to successfully characterize
various composite materials. Derivation of the model parameters from the experiments is not always
straightforward. Derivation of Material Properties (p. 43) outlines the required steps in extracting the
parameters from experiment into a form suitable for input into Autodyn.

Example applications are presented in Example Applications (p. 53) and recommendations on how to
perform a typical composite analysis are discussed in Recommendations (p. 63).

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 1
 Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
2 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 2: Principal Directions in Autodyn
An orthotropic material has properties that are different in three mutually perpendicular directions and
has three mutually perpendicular planes of material symmetry. Thus the properties of such a material
are a function of orientation of the material with respect to the global co-ordinate system. The ortho-
tropic constitutive models to be described in 3 are done so in terms of the principal material directions.
It is therefore required to define the initial orientation of these principal material directions in an
Autodyn model with respect to the global co-ordinates. Definition of the principal material direction is
described in Autodyn-2D (p. 3) for Autodyn-2D and in Autodyn-3D Lagrange and ALE Parts (p. 5)
and Autodyn-3D Shell Parts (p. 7) for Autodyn-3D.

2.1. Autodyn-2D
For Lagrange and ALE subgrids, two of the principal directions lie in the XY plane. For planar symmetry
the third principal direction is perpendicular to the XY plane, and for axial symmetry in Autodyn-2D
the third principal direction is the hoop direction (this direction is termed as either TT or 33 in Autodyn-
2D).

The first principal direction is defined as the angle, in degrees, between the principal direction and the
global XY axes. This angle, which can be examined or plotted as a time-history, is called "[Link]".
The second principal direction is orthogonal to the first principal direction in the XY plane. You can plot
the principal directions using the "Direction" menu option.

You are allowed to define the initial orientation of the principal material axes in one of three ways. In
all three cases, a rotation angle, β, can be specified. If this angle is non-zero, the first and second prin-
cipal axes are additionally rotated anti-clockwise through this angle.

2.1.1. X-Y Space

The first principal axis is defined relative to X-axis.

Figure 2.1: Material directions in X-Y space

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 3
Principal Directions in Autodyn

2.1.2. Polar Space

The first principal axis for each cell is relative to the direction of the line drawn from the polar origin
to the center of the cell.

Figure 2.2: Material directions in polar space

2.1.3. I-J-K Space

The first principal axis for each cell is relative to the direction of increasing I for that cell.

Figure 2.3: Material directions in I-J space

where,

(2.1)

Note:

When the I-J-K Space option is selected for 2D unstructured Lagrange elements, the direction
of the first principal axis will coincide with the X-axis direction.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
4 of ANSYS, Inc. and its subsidiaries and affiliates.
 Autodyn-3D Lagrange and ALE Parts

2.2. Autodyn-3D Lagrange and ALE Parts

For 3D Lagrange and ALE subgrids, the three principal directions are oriented in XYZ space and are
mutually perpendicular. These three directions are defined by the use of four variables as explained
below, and illustrated in Figure 2.4: Material directions in Autodyn-3D defined in XYZ-space (p. 5).

The first principal direction is a vector in XYZ space (1 in Figure 2.4: Material directions in Autodyn-3D
defined in XYZ-space (p. 5)) and is defined in Autodyn-3D by the x, y and z components of this vector.
The second and third principal directions are two vectors (2 and 3) which lie in a plane that has it’s
normal as the first principal direction. The location of the second and third principal directions are
defined in Autodyn-3D by a single angle; this angle (α) is defined as being that between the second
principal direction and the vector (2’) which is the cross-product (in other words, is normal to) the X
axis and the first principal direction. In Figure 2.4: Material directions in Autodyn-3D defined in XYZ-
space (p. 5), as an example, the direction 1 is shown as vector which is rotated in the XY plane by an
angle θ; therefore the direction 2’ is coincident with direction Z. The direction 2 is rotated from 2’ by
an angle α in the 23 plane. Note that the principal angle can vary between -90° and +90°. Direction 3
is normal to directions 1 and 2.

Figure 2.4: Material directions in Autodyn-3D defined in XYZ-space

You are allowed to define the initial orientation of the principal material axes in one of two ways, as
shown below:

2.2.1. X-Y-Z Space

The first principal direction is defined relative to X-Y-Z axes using the x, y and z components of the
vector defining the direction. Directions 2 and 3 are defined using the angle α as explained above.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 5
Principal Directions in Autodyn

Figure 2.5: Material directions in Autodyn-3D defined in XYZ-space

2.2.2. I-J-K Space

The first principal direction for each cell is relative to the direction of increasing I for each cell, as
shown below. The directions are calculated individually for each cell when the material is filled into
the cells; then the directions are stored as four variables exactly as for the X-Y-Z space option. The
directions are found as follows:

• The direction 1 is a vector from the centroid of the cell face I-1 (point a) to the centroid of cell face I (point
b)

• Next find the midpoints of the K-lines at [I-1,J-1] and [I,J], which are points p and q respectively

• Direction 3 is the cross product of direction 1 and the vector through points p and q

• Direction 2 is the cross product of directions 1 and 3

Using the I-J-K option you can define initial material directions that are aligned with the cells. This
can be useful when modelling shapes such as cylinders where the principal directions are aligned
with the major axes of the body.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
6 of ANSYS, Inc. and its subsidiaries and affiliates.
 Autodyn-3D Shell Parts

Figure 2.6: Material directions in Autodyn-3D defined in IJK-space

Note:

When the I-J-K Space option is selected for 3D unstructured Lagrange elements, the direction
of the first principal axis will coincide with the global X-axis direction.

2.3. Autodyn-3D Shell Parts

For 3D shell parts, a local co-rotational element coordinate system is used for the geometrical definition
of the shell element.

The local co-rotational coordinate directions for structured shells are found as follows:

• The direction 1 is a vector from the centroid of the element face K (point a) to the centroid of the element
face K-1 (point b).

• Next find the centroids of the faces at J and J-1, which are points p and q.

• Direction 3 is the cross product of direction 1 and the vector through points p and q.

• Direction 2 is the cross product of directions 1 and 3.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 7
Principal Directions in Autodyn

Figure 2.7: Co-rotational coordinate system used for 3D structured Shell elements

The local co-rotational coordinate directions for unstructured shells are found in a similar fashion as for
structured shells and are shown in Figure 2.8: Co-rotational coordinate system used for 3D unstructured
Shell elements (p. 8):

Figure 2.8: Co-rotational coordinate system used for 3D unstructured Shell elements

The first principal material direction is defined by a vector in XYZ space (XYZ in Figure 2.7: Co-rotational
coordinate system used for 3D structured Shell elements (p. 8) and Figure 2.8: Co-rotational coordinate
system used for 3D unstructured Shell elements (p. 8)) and makes an angle (  ) with the first co-rota-
tional direction 1.

For a 3D shell element, the first material direction (11) lies in the 1-2 co-rotational coordinate plane and
its direction is obtained through a rotation from the co-rotational 1 direction over an angle    around
the co-rotational 3 direction (Figure 2.9: Material directions defined in XYZ-space for 3D Shell ele-
ments (p. 9)).

The third material direction 33 coincides with the co-rotational direction 3. The second material direction
(22) is defined in the 1-2 co-rotational coordinate plane, and is normal to the plane through the first
material direction 11 and the co-rotational 3 direction.

Note:

Angle    can vary between 0° and 180°.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
8 of ANSYS, Inc. and its subsidiaries and affiliates.
 Autodyn-3D Shell Parts

Figure 2.9: Material directions defined in XYZ-space for 3D Shell elements

You are allowed to define the initial orientation of the principal material axes in one of two ways, as
shown in X-Y-Z Space (p. 9) and I-J-K Space (p. 9).

2.3.1. X-Y-Z Space

The first principal direction is defined relative to X-Y-Z axes using the x, y and z components of the
vector defining the direction. Material directions 11 and 22 are defined using the angle    as explained
in Figure 2.9: Material directions defined in XYZ-space for 3D Shell elements (p. 9).

2.3.2. I-J-K Space

The first principal direction for each structured shell element is relative to the direction of decreasing
K for each element, as shown in Figure 2.10: Material directions defined in IJK-space for 3D structured
Shell (p. 9). The material directions will coincide with the co-rotational directions of the shell element.

The material directions will also coincide with the co-rotational directions of the shell element when
the I-J-K option is used for unstructured shell elements.

Using the I-J-K option you can define initial material directions that are aligned with the shell elements.
This can be useful when modelling shapes such as cylinders where the principal directions are aligned
with the major axes of the body.

Figure 2.10: Material directions defined in IJK-space for 3D structured Shell

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 9
 Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
10 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 3: Orthotropic Constitutive Models
In general the behavior of composite laminates can be represented through a set of orthotropic con-
stitutive relations. In Orthotropic Elastic Model (p. 11) a set of such relationships are described which
assume the material behavior remains elastic and the volumetric response linear. For more complicated
material response a methodology was developed, under contract from ESA (European Space Agency)
[1], which allows a nonlinear equation of state to be used in conjunction with an orthotropic stiffness
matrix and is described in Equations of State (p. 14). This is important when modeling applications such
as hypervelocity impacts.

Some composite materials, such as Kevlar-epoxy, exhibit significant nonlinear stress-strain relationships
(see Figure 3.1: Typical In-Plane Stress-Strain Behavior of Kevlar-Epoxy (p. 11)). In order to model such
observed nonlinear behavior an orthotropic hardening model has been implemented. Again developed
under contract from ESA [2] it uses an anisotropic plasticity based loading/failure surface and is described
in Orthotropic Yield Strength and Hardening Model (p. 17).

Figure 3.1: Typical In-Plane Stress-Strain Behavior of Kevlar-Epoxy

3.1. Orthotropic Elastic Model

For conventional materials the constitutive relations, or strength model, defines the deviatoric response
of the material as a function of any combination of deviatoric strain, deviatoric strain rate, temperature,
pressure, and damage.
(3.1)

For many other materials, including composites, the macroscopic properties are not identical in all dir-
ections. In general, the behavior of such materials is represented through a set of orthotropic constitutive
relations. Constitutive relations for this type of material are conventionally based on a total stress for-
mulation, as opposed to dividing the total stress into hydrostatic and deviatoric components. Thus, the
incremental stress-strain relations can be expressed as

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 11
Orthotropic Constitutive Models

(3.2)

where

[C] = stiffness matrix

= strain rate tensor

= time step

The linear elastic constitutive relations for a general anisotropic (triclinic) material can be expressed in
relation to a Cartesian co-ordinate system, in contracted notation, as

(3.3)

In which there are 21 independent elastic constants, Cij . If there is one plane of material symmetry, the
above stress strain relations reduce to

(3.4)

where the plane of symmetry is x33 = 0. Such a material is termed monoclinic and there are 13 inde-
pendent elastic constants. This is the basic form assumed for textiles and composites which generally
have a plain of symmetry in the through thickness direction while the symmetry in the plane of the
material is dependant on the lay-up and weave of the fibers.

If a second plane of symmetry exists in the plane of the fiber-composite then symmetry will also exist
in a third mutually orthogonal plane. Such a material is said to be orthotropic and the constitutive rela-
tions are of the form:

(3.5)

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
12 of ANSYS, Inc. and its subsidiaries and affiliates.
 Orthotropic Elastic Model

The inverse of the above stiffness matrix for a three-dimensional orthotropic configuration, the compliance
matrix, is

(3.6)

where,

Eii are the Young’s Moduli in the principal material directions,

Gij are the shear moduli,

are the Poisson’s ratios, where is defined as the transverse strain in the j-direction when stressed
in the i-direction, that is:

; note that

Note:

In Autodyn 2D the directions 23 and 31 do not apply, so:

G23 and G31 are not required

= =0

σ23= σ31= 0

The above, and similar, constitutive relations for orthotropic materials are commonly applied to fiber
reinforced composites in which the fibers are set in a solid matrix material. Inherent in the model is the
assumption of a linear volumetric elastic response of the material. This may not be representative of
actual material response under the high pressures experienced during a hypervelocity impact event.

There are, of course, restrictions on the elastic constants that can be used for an orthotropic material
and these restrictions are more complex than those for isotropic materials. These restrictions result from
the fact that the sum of work done by all stress components must be positive in order to avoid the
creation of energy. The first condition states that the elastic constants are positive:
E11 ,E22 ,E33 ,G12 ,G23 ,G31 > 0 (3.7)

Second, the determinant of the stiffness matrix should be positive:

(3.8)

Finally the requirement for positive stiffnesses leads to:

, , (3.9)

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 13
Orthotropic Constitutive Models

Whenever the elastic constants in an Autodyn model are defined or redefined the three conditions
above are tested and the user is informed if any of them are violated.

3.2. Equations of State

The material model developed in [1] provides a mechanism in an orthotropic material to calculate:

• The contributions to pressure from the isotropic and deviatoric strain components

• The contributions to the deviatoric stress from the deviatoric strains

Further, this methodology gives rise to the possibility for incorporating nonlinear effects (such as shock
effects) that can be attributed to the volumetric straining in the material. To use this model ‘Ortho’ is
selected as the equation state for the material and either ‘Polynomial’ or ‘Shock’ for the volumetric re-
sponse option.

The incremental linear elastic constitutive relations for an orthotropic material can be expressed, in
contracted notation, as:

(3.10)

In order to include nonlinear shock effects in the above linear relations, it is first desirable to separate
the volumetric (thermodynamic) response of the material from its ability to carry shear loads (strength).
To this end, it is convenient to split the strain increments into their average, , and deviatoric, ,
components.
(3.11)

Now, defining the average direct strain increment, , as a third of the trace of the strain tensor,
(3.12)

and assuming, for small strain increments, the volumetric strain increment is defined as
(3.13)

The total strain increments can be expressed in terms of the volumetric and deviatoric strain increments
resulting in the following orthotropic constitutive relation.

(3.14)

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
14 of ANSYS, Inc. and its subsidiaries and affiliates.
 Equations of State

If the above relations are expanded and the deviatoric and volumetric terms grouped, the following
expressions for the direct stress increments results.

(3.15)

To find the equivalent pressure increment, we first define the pressure as a third of the trace of the
stress increment tensor;
(3.16)

Substituting Equation 3.15 (p. 15) into Equation 3.16 (p. 15) results in an expression for the pressure
increment of the form

(3.17)

from which the contributions to the pressure from volumetric and deviatoric components of strain can
clearly be identified.

For an isotropic material, the stiffness matrix coefficients can be represented in terms of the material
bulk modulus, K, and shear modulus, G. Thus,

(3.18)

Substituting Equation 3.18 (p. 15) into Equation 3.16 (p. 15) gives
(3.19)

and given
(3.20)

Equation 3.19 (p. 15) reduces to

(3.21)

which is immediately recognizable as the standard relationship between pressure and volumetric strain
(Hooke’s law) at low compressions.

The first term of Equation 3.17 (p. 15) can therefore be used to define the volumetric (thermodynamic)
response of an orthotropic material in which the effective bulk modulus of the material K’ is
(3.22)

For the inclusion of nonlinear shock effects, the contribution to pressure from volumetric strain is
modified to include nonlinear terms. The final incremental pressure calculation becomes

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 15
Orthotropic Constitutive Models

(3.23)

where the pressure contribution ΔPEOS from volumetric strains can include the nonlinear shock (ther-
modynamic) effects and energy dependence as in a conventional equation of state.

A form of equation of state that is used extensively for isotropic solid continua is known as the Mie-
Grüneisen form:
(3.24)

Where the Grüneisen gamma is defined as

(3.25)

The functions pr(v) and er(v) are assumed to be known functions of v on some reference curve.

Two Mie-Grüneisen forms of equation of state are available for coupling with an orthotropic response
in the AMMHIS model and are now described.

3.2.1. Shock Equation of State

To define the pressure and energy on the Hugoniot this equation of state requires the material shock
velocity-particle velocity relationship as material input. In Autodyn this is specified through
Us = C0 + Sup (3.26)

where

Us is the shock velocity

up is the particle velocity
S is the slope of the Us – up relationship
C0 is the bulk acoustic sound speed

For the case of an orthotropic material, the bulk acoustic sound speed is calculated from the effective
bulk modulus Equation 3.22 (p. 15) and density:
(3.27)

3.2.2. Polynomial Equation of State

This is of the form
(3.28)

where the first term in (3-28) is equivalent to a linear equation of state with the bulk modulus, K’,
derived from the orthotropic material stiffness coefficients.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
16 of ANSYS, Inc. and its subsidiaries and affiliates.
 Strength Models

3.3. Strength Models

Three strength models are available for orthotropic materials in Autodyn-2D and -3D, as follows:
3.3.1. Elastic
3.3.2. VonMises
3.3.3. Orthotropic Yield Strength and Hardening Model

3.3.1. Elastic
The orthotropic material will respond elastically to loading until failure is reached.

3.3.2. VonMises
The orthotropic material initially will respond elastically to loading until a stress state is reached that
exceeds the constant VonMises yield value and the material starts to deform plasticly until failure is
reached.

3.3.3. Orthotropic Yield Strength and Hardening Model

In this section we describe in detail the implementation of a quadratic limit surface used to model
the anisotropic hardening behavior observed in some composite materials.

There are perhaps two definitions for failure; for a structural designer failure may be considered the
point at which the material starts to become nonlinear (due to plasticity or micro-cracking, for example),
however for the simulation of extreme loading events such as HVI failure is considered to be the
point at which the material actually ruptures (becomes perforated). To get to this point, the material
transitions elastic to inelastic deformation (due to micro-cracking, plasticity in the matrix or re-orient-
ation of fibers) and finally reaches ultimate failure. We take the later definition, and for the purposes
of this model we will term this nonlinear stress strain behavior “hardening” (see Figure 3.1: Typical
In-Plane Stress-Strain Behavior of Kevlar-Epoxy (p. 11)).

The anisotropic yield criteria of Tsai-Hill can be used to represent the potentially different limits in
elastic material behavior in the three orthotropic material directions [5]. This type of surface is often
used to model anisotropic yielding in materials with isotropic stiffness such as metals and polymers.
This method does however have some limitations for materials with orthotropic stiffness. The first is
that it is assumed that all inelastic deformation is at constant volume; in other words, due to deviat-
oric strain. However, for anisotropic materials, deviatoric strain contributes to the volumetric pressure.
Hence it may be inconsistent to not change the volume (and therefore pressure) when inelastic de-
formation takes place. The model also does not include any hardening, the yield/initial failure surface
remains at a constant level of stress.

A nine-parameter yield function was presented by Chen et al [3] and has been implemented into the
Autodyn code [2]. This model relaxes the constant pressure assumption in Hill’s theory and therefore
is more generally applicable to materials with orthotropic stiffness. With careful selection of the
parameters one can however return to the constant volume assumption. The yield function also in-
cludes a hardening parameter. Softening behavior and damage is not included with this yield function
but it should be generally applicable to fiber reinforced composite materials. Also, because of the
generality of the yield function it is applicable to both isotropic and anisotropic materials. The addi-
tional benefit of using a plasticity based hardening approach is that it can be used in conjunction
with the existing nonlinear shock response features described in Equations of State (p. 14).

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 17
Orthotropic Constitutive Models

[Link]. Quadratic Limit Surface

The quadratic yield/flow surface [3] was selected to represent nonlinear hardening effects.

(3.29)

The yield function is quadratic in material stress space and includes nine material constants, , to
represent the degree of anisotropy in the material behavior. The parameter k varies with the effective
inelastic strain in the material and can be used to represent hardening behavior.

This yield surface is very general and by careful selection of the coefficients can reduce to several
other well-known yield criteria. The yield criteria Equation 3.29 (p. 18) reduces to Hill's orthotropic
yield function under the following conditions:

(3.30)

Also, the von Mises J2 yield criteria is recovered if the following values are set for the plasticity
parameters:

(3.31)

The plasticity parameters would ideally be calibrated from the experimental stress-strain data ob-
tained from three simple uniaxial tension tests and three pure shear tests. The six parameters asso-
ciated with the normal stresses are then determined from the definition of a master effective stress-
effective plastic strain relationship and from the definition of Plastic Poisson’s’ ratios (PPR).
The three constants associated with the shear stresses are then obtained by collapsing the corres-
ponding data onto the curve. In practice all the necessary stress-strain data are not usually
available and assumptions have to be made. Indeed in [3] the required stress-strain data was obtained
by 3D micro mechanical simulations.

It is necessary that the plasticity parameters, to , define a real closed surface in stress space.
To ensure this requirement is met the following constraints are placed on the plasticity parameters,

>0

Additionally, the following two matrices are defined,

where it is required that Det E < 0, and, the non-zero eigen values of the matrix e all have the same
sign.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
18 of ANSYS, Inc. and its subsidiaries and affiliates.
 Strength Models

After initial yielding material behavior will be partly elastic and partly plastic. In order to derive the
relationship between the plastic strain increment and the stress increment it is necessary to make
a further assumption about the material behavior. The incremental plastic strains are defined as
follows,

(3.32)

These are the Prandtl-Reuss equations, often called an associated flow rule, and state that the plastic
strain increments are proportional to the stress gradient of the yield function. The proportionality
constant, dλ , is known as the plastic strain-rate multiplier. Written out explicitly the plastic strain
increments are given by

(3.33)

PPR, , are defined as,

, (no summation over repeated indices) (3.34)

From Equation 3.33 (p. 19) and Equation 3.34 (p. 19) the following relationships are easily derived

(3.35)

The incremental effective plastic strain can be calculated through a concept of plastic work,
(3.36)

We can rewrite Equation 3.36 (p. 19) as,

(3.37)

If the plastic strain increment is replaced by the flow rule Equation 3.32 (p. 19) and multiplied by
the transpose of the stress tensor, and using the following definition of effective stress,
(3.38)

it is easily shown that the effective plastic strain increment for the quadratic yield function is,
(3.39)

The following explicit definition of incremental effective plastic strain results from substitution of
Equation 3.29 (p. 18) into Equation 3.39 (p. 19) and using Equation 3.33 (p. 19)

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 19
Orthotropic Constitutive Models

(3.40)

where

and

[Link]. Stress Return Algorithm

The typical series of calculations that are carried out each incremental time step (or cycle) in an
element are shown schematically in Figure 3.2: Computational Cycle (p. 20). Starting at the bottom
of Figure 3.2: Computational Cycle (p. 20) the boundary and/or interactive forces are updated and
combined with the forces for inner zones computed during the previous time cycle. Then, for all
non-interactive nodes, the accelerations, velocities and positions are computed from the momentum
equation and a further integration. From these values the new zone volumes and strain rates may
be calculated. With the use of a material model together with the energy equation the zone pres-
sures, stresses and energies may be calculated, thus providing forces for use at the start of the next
integration cycle.

Figure 3.2: Computational Cycle

Having determined the strain rates and the volume change the stresses can be calculated.
These are then checked against the quadratic limit criteria of equation Equation 3.29 (p. 18). Initially

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
20 of ANSYS, Inc. and its subsidiaries and affiliates.
 Strength Models

all stresses are updated using an elastic relationship. If these stresses remain within the yield surface
then the material is assumed to have either loaded or unloaded elastically. If, however, the resultant
updated stress state lies outside of the yield surface, point B in Figure 3.3: Yield surface (p. 21),
steps must be taken to return the stresses normal to the yield surface.

Figure 3.3: Yield surface

The backward-Euler stresses are calculated as,

(3.41)

where, C is the elastic stiffness matrix and from the flow rule,

(3.42)

which involves the vector aC that is normal to the yield surface at the final position C. At this location
the stress state satisfies the yield criteria. Except in special cases this vector aC cannot be determined
from the data at point B. Hence an iterative procedure must be used.

As a first step a predictor stress state, is calculated using the following equation,
(3.43)

where,

(3.44)

dλ here is defined so as to make the value of the yield function zero when a Taylor expansion of
the yield criteria is performed about the trial point B. In general this will result in a stress state that
still remains outside of the yield surface. This is because the normal at trial elastic state B is not
equal to the final normal at point D (see Figure 3.4: Schematic representation of backward-Euler
return algorithm (p. 22)). Further iterations are therefore usually required.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 21
Orthotropic Constitutive Models

Figure 3.4: Schematic representation of backward-Euler return algorithm

The stress return is achieved using a backward-Euler return algorithm and the yield function is as-
sumed to be satisfied if the returned stress state is within 1% of the yield surface.

Once the yield function is satisfied the plastic portion of the strain increment is determined from
the following,
(3.45)

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
22 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 4: Orthotropic Material Failure Models
The following orthotropic material failure models are described in this chapter:
4.1. Brittle Damage Model
4.2. Orthotropic Damage Model

4.1. Brittle Damage Model

For brittle materials, three orthotropic failure initiation models are available in Autodyn-2D and 3D, as
follows:

• Material Stress

• Material Strain

• Material Stress/Strain.

These models allow different tensile and shear failure stresses and/or strains for each of the principal
(material) directions. The models are intended to enable modeling of orthotropic failure.

Importantly, the orthotropic failure models can be used together with the orthotropic equation of state,
which before failure uses the incremental elastic stress-strain relations detailed in Orthotropic Constitutive
Models (p. 11). With any of the above Material failure models, failure is initiated in a principal direction
when the stress or strain reaches a user-specified limiting value.

4.1.1. Failure Initiation

Material Stress
This model is useful for materials that are likely to fail along predefined material planes; for instance,
where delamination failure occurs between layers in a composite plate.

Tensile Failure Stress 11

Tensile Failure Stress 22
Tensile Failure Stress 33
Maximum Shear Stress 12
Maximum Shear Stress 23 (3D only)
Maximum Shear Stress 31 (3D only)

Failure is initiated if any of the principal material stresses exceed their respective tensile failure
stresses. For shear failures the shear stress on planes parallel to the principal directions are checked
against the maximum shear stress.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 23
Orthotropic Material Failure Models

Material Strain
This model is useful for materials, which are likely to fail along predefined material planes: For instance,
where failure occurs parallel to the interface between two layers in a laminated plate.

Tensile Failure Strain 11

Tensile Failure Strain 22
Tensile Failure Strain 33
Maximum Shear Strain 12
Maximum Shear Strain 23 (3D only)
Maximum Shear Strain 31 (3D only)

Failure is initiated if any of the principal material strains exceed their respective tensile failure
strains. For shear failures the shear strain on planes parallel to the principal directions are checked
against the maximum shear strain.

Material Stress/Strain
This model is useful for materials that are likely to fail along predefined material planes: For instance,
where failure occurs, parallel to the interface, between two layers in a laminated plate.

Tensile Failure Stress 11

Tensile Failure Stress 22
Tensile Failure Stress 33
Maximum Shear Stress 12
Maximum Shear Stress 23 (3D only)
Maximum Shear Stress 31 (3D only)
Tensile Failure Strain 11
Tensile Failure Strain 22
Tensile Failure Strain 33
Maximum Shear Strain 12
Maximum Shear Strain 23 (3D only)
Maximum Shear Strain 31 (3D only)

Failure is initiated if any of the principal material stresses or strains exceed their respective failure
levels. For shear failures the shear stresses, and strains, on planes parallel to the principal directions
are checked against the inputted maximum shear stresses, and strains.

4.1.2. Post-Failure Response

When an orthotropic equation of state is used in conjunction with material stress, strain or stress/strain
failure criteria an option to use isotropic or orthotropic post failure response is included in the failure
data input. Subsequent to failure initiation, the failed cell stiffness and strength properties are modified
depending on the failure initiation modes described in the subsections below.

[Link]. Isotropic Post-Failure

The failure models outlined in Failure Initiation (p. 23) allow different tensile and shear failure
stresses and/or strains for each of the “principal” directions. If the isotropic post-failure option is

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
24 of ANSYS, Inc. and its subsidiaries and affiliates.
 Brittle Damage Model

selected failed cells can only carry bulk compressive stresses. Therefore after failure is initiated in
a cell the following occurs:

• The principal (material) stress in the direction of failure is set to zero.

• All shear moduli are set to zero.

• All shear stresses are set to zero.

• The average stress (in other words, pressure) is recomputed, using the normal calculation:

For the orthotropic equation of state, post-failure behavior is modelled as detailed below. This post-
failure model is in effect an isotropic post-failure response.

• The orthotropic elastic incremental stress-strain relations are applied.

• The average stress (in other words, pressure) is recomputed, using the calculation above.

• If the cell is in compression,

The principal stresses are set equal to the average stress (in other words, pressure):

• If the cell is in tension,

All principal stresses, and therefore the average stress (in other words, pressure), are set to zero:

[Link]. Orthotropic Post-Failure Response

The post-failure behavior of laminated composite materials is better represented through the ortho-
tropic post failure option. This option was developed specifically for simulating the performance
and failure, including delamination, of fiber reinforced composite materials.

If the orthotropic option is selected, the user is prompted to define the post-failure parameters.
This includes the post failure response mode for failure in each material direction and the failed
material residual shear strength (see Table 4.1: Orthotropic Post-Failure Options (p. 25)).

Table 4.1: Orthotropic Post-Failure Options

Failed Direction Post Failure Post Failure Response

Option
FAIL 11 11 only Zero tensile stress in 11 direction.

Bulk Zero tensile stress in all directions.

FAIL 22 22 only Zero tensile stress in 22 direction.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 25
Orthotropic Material Failure Models

Failed Direction Post Failure Post Failure Response

Option
Bulk Zero tensile stress in all directions.
FAIL 33 33 only Zero tensile stress in 33 direction.

Bulk Zero tensile stress in all directions.

FAIL 12 12 & 11 Zero tensile stress in 11 direction.
only
Zero tensile stress in 22 direction.
12 & 22
only Zero tensile stress in all directions.

Bulk
Residual Shear 0.0 to 1.0 Residual shear modulus is set to the specified
Stiffness Fracture value times the intact shear modulus. Default
value is 0.2.
Maximum 0.0 to 1.0e20 Maximum shear stress allowed in a failed cell. A
Residual Shear value equal to, or less than, the failure shear
Stress stress is recommended.

For all failure modes, on failure initiation, the stress in the failed material directions are set to zero.
In addition the stresses in material directions orthogonal to the failed direction are reduced to ac-
count for the loss in the Poisson effect from strain in the failed direction. In subsequent cycles, cell
tensile stresses are only allowed in non-failed directions.

If a cell fails in two or more directions, bulk failure is assumed.

Subsequent to failure initiation, the failed cell stiffness and strength properties are modified depend-
ing on the failure initiation modes described below.

Delamination
In an axisymmetric simulation, the 11-direction is assumed to be through the thickness of the laminate
and the 33-direction is the hoop direction. Delamination can result from excessive through thickness
tensile stresses and/or strains or from excessive shear stress and or strain in the 12 plane. If failure is
initiated in either of these two modes, the stress in the 11-direction is instantaneously set to zero and
the strain in the 11-direction at failure is stored. Subsequently, if the tensile material strain in the 11-
direction exceeds the failure strain, the material stiffness matrix is modified as

(4.1)

This stiffness modification does not allow tensile through thickness stresses while tensile in-
plane stresses are maintained.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
26 of ANSYS, Inc. and its subsidiaries and affiliates.
 Brittle Damage Model

Also note that the delamination will in practice be associated with a reduction in shear stiffness.
Often, in the absence of appropriate material data, a nominal value of 20% is typically used for
the residual shear stiffness α.

In-plane Failure
In an axisymmetric simulation, the 22- and 33-directions are assumed to be in the plane of the com-
posite (in other words, in the fiber directions). If failure is initiated in these two modes, the stress in
the failed direction is instantaneously set to zero and the strain in the failed direction at failure is
stored. Subsequently, if the tensile material strain in the failed direction exceeds the failure strain, the
material stiffness matrix is modified as

22-failure:

(4.2)

33-failure:

(4.3)

This stiffness modification does not allow tensile stresses in the failed directions. Also note that
these failure modes will in practice be associated with a reduction in shear stiffness. Often, in
the absence of appropriate material data, a nominal value of 20% is typically used for the residual
shear stiffness α.

Combined Failure
The combined effect of failure in all three material directions is represented by a change in the mater-
ial stiffness and strength to isotropic with no stress deviators and no tensile material stresses.

Melting and Decomposition

Melting/vaporisation of epoxy and decomposition of the fiber material has been observed in composites
with an epoxy matrix, for example Kevlar-epoxy [1]. This occurred in a finite region directly under the
impact point under uniaxial strain conditions and velocities above or around 1000m/s.

To represent this phenomena in an approximate way in a numerical model, the following features
are available. An epoxy melting temperature can be specified. It is assumed that this has a very
similar effect to delamination in the laminate. The procedure outlined in Section [Link].1 is
followed for this failure initiation mode.

Additionally, a decomposition temperature for the fiber can be specified. Decomposition may
lead to an inhomogeneous material of unknown properties. The model therefore assumes that
the decomposed material has properties of the intact material under bulk compression. In bulk
tension, pressure, deviatoric and tensile stresses are set to zero.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 27
Orthotropic Material Failure Models

4.2. Orthotropic Damage Model

Composite materials are not perfectly brittle, instantaneous failure is not necessarily observed in all or
any of the material directions. The stress/strain based failure models available for use in conjunction
with the orthotropic model instantaneously set material stresses to zero on failure initiation, as described
in Brittle Damage Model (p. 23).

In this section we report developments [2] which model the softening behavior observed in some
composites (see Figure 3.1: Typical In-Plane Stress-Strain Behavior of Kevlar-Epoxy (p. 11)). Maximum
tensile and shear stresses in a cell are limited by failure surfaces used to define failure initiation. Different
failure modes are considered each of which are described by a unique surface. To model the gradual
reduction in the ability of a laminate to carry stress a crack softening approach is used. This reduces
the maximum tensile stress that can be sustained in an element as a function of some measure of crack
strain.

4.2.1. Failure Initiation

Hashin failure criteria have been extensively used in modelling the prediction of impact damage in
composites. However, in these criteria for fiber and matrix failure only plane stresses, σ22, σ23 and
σ33, are considered. As a consequence the applicability of these criteria for out of plane impacts is
questionable.

Modified versions of these failure criteria along with a criterion for delamination are presented in [6]
and have been implemented into Autodyn [2]. For the fiber failure and matrix cracking criteria out
of plane shear stresses are included in addition to the original criteria. The failure initiation criteria
are listed below.

Delamination
(4.4)

Fiber Failure
(4.5)

Matrix Cracking
(4.6)

The initiation criteria Equation 4.4 (p. 28) to Equation 4.6 (p. 28) are only applied in tension.

4.2.2. Damage Model

The failure initiation criteria can be considered as surfaces that define regions of valid stress. If a state
of stress is found to lay outside of this surface a backward-Euler algorithm is used to return the stress
to the failure surface. The resulting inelastic increment in strain is then accumulated as crack strain
εcr . The maximum stress that can be sustained in an element is then reduced as a function of crack
strain.

In this approach the area under the softening portion of the stress/strain curve is related to the fracture
energy Gf , which is a material property. This is simply illustrated in Figure 4.1: Schematic Illustration
of Crack Softening Algorithm (p. 29).

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
28 of ANSYS, Inc. and its subsidiaries and affiliates.
 Orthotropic Damage Model

Figure 4.1: Schematic Illustration of Crack Softening Algorithm

A mathematical description of Figure 4.1: Schematic Illustration of Crack Softening Algorithm (p. 29)
can be written as,

(4.7)

where L is a characteristic cell dimension in the direction of failure and is included to improve the
objectivity of the solution.

Failure is initiated when the stress reaches the value required for failure σfail . At this point the crack
strain εcr is zero. A linear softening slope is assumed and therefore the ultimate crack strain εU , the
strain at which tensile stresses can no longer be sustained, is calculated as:

(4.8)

The gradient of the linear softening slope is given by:

(4.9)

After failure initiation a linear damage law is used to reduce the maximum material stress in a cell as
a function of crack strain.
(4.10)

At initiation εcr = 0, therefore Dam = 0 . When the cell has no strength, Dam = 1. At any point between
these times the maximum tensile material stress that can be supported in the cell is:
(4.11)

Rather than define an effective crack strain analogous to effective plastic strain, damage has been
implemented as a full tensor

(4.12)

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 29
Orthotropic Material Failure Models

Each component of Equation 4.12 (p. 29) has an independent damage variable associated with it
resulting in orthotropic damage. If we examine the 11-plane softening surface, damage is included
for each stress component as follows,

(4.13)

As with the hardening behavior of Orthotropic Yield Strength and Hardening Model (p. 17) an asso-
ciative flow rule is assumed when returning the stresses to the softening surfaces,
(4.14)

Considering the flow rule Equation 4.14 (p. 30) the incremental crack strain components are calculated
as follows for the 11-plane softening surface Equation 4.13 (p. 30),

(4.15)

The possibility that an increase in damage in a particular direction may reduce the maximum strength
in other directions can be incorporated via the damage coupling coefficient, C. The value of C can
vary between 0 and 1. For example, for the failure surface considered above, the damage increments
are updated as,

Similar equations are trivially derived for the other softening surfaces.

If an element is fully damaged in two or more material directions, bulk failure is assumed:

Damage Post Failure Response

Direction
D11 = 1.0
D22 = 1.0 Zero tensile stress in all directions

D22 = 1.0
D33 = 1.0 Zero tensile stress in all directions

D11 = 1.0
Zero tensile stress in all directions
D33 = 1.0

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
30 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 5: Material Characterization Tests
To fully characterize a composite material requires an extensive series of experimental tests. The number
of tests required though may vary depending on the material behavior and the application. This section
briefly describes experiments that have been used to obtain composite material properties for studies
using the material models described in 3 and 4. The testing strategy described in this chapter has been
developed under ESA research contract No. 12400/97/NL/PA(SC) and performed by Ernst Mach Institute
(EMI), Freiburg, in conjunction with Century Dynamics, [1][2].

5.1. Directional Strength Properties

Normal in-plane testing of composites in compression and tension is standardized: for example standard
EN ISO 527-4. Tensile and shear tests measure important material characteristics such as elastic and
plastic Poisson ratios, nonlinear hardening properties and failure limits, and are discussed in the following
sections.
5.1.1. In-Plane Tension Tests
5.1.2. Out-of-Plane Shear

5.1.1. In-Plane Tension Tests

0° Tension Test
Tensile tests are typically performed to determine the stress-strain response of a composite material in
the plane of the laminate. Data commonly generated by this test are the stiffness and strength properties.
Samples, of similar shape to that illustrated in Figure 5.1: Schematic Diagram of a Tensile Test Speci-
men (p. 32), are cut from the laminate sheets and the ends of the sample are clamped and separated
at a constant speed in opposite directions using a machine such as an Instron®. The tension or force re-
quired to do this is recorded. The test results are shown as a stress-strain curve, the shape of which
provides information about the material behavior of the sample.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 31
Material Characterization Tests

Figure 5.1: Schematic Diagram of a Tensile Test Specimen

Strains can be recorded by using the crosshead displacement of the experimental apparatus, or
using an extensometer. Alternatively biaxial strain gauges located both on the side and front of
the test specimen allow for measurement of the strains in the principal material directions. In this
way both in-plane and through thickness properties can be measured. For unidirectional laminates,
samples can be prepared such that the specimens can be orientated either in the direction of the
fibers (0° tension test) or at 90° to the fibers (90° tension test).

Figure 5.2: Calculation of In Plane and Out of Plane Poisson Ratios (p. 32) and Figure 5.3: Typical
Recorded Uniaxial Stress-Strain Relationship (p. 33) show a typical result of a tensile test on a
composite material.

Figure 5.2: Calculation of In Plane and Out of Plane Poisson Ratios

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
32 of ANSYS, Inc. and its subsidiaries and affiliates.
 Directional Strength Properties

Figure 5.3: Typical Recorded Uniaxial Stress-Strain Relationship

45° Tension Test

Tension tests on samples manufactured so that the fibers are aligned at 45° orientation to the applied
load are performed to derive the in-plane shear modulus according to Equation 5.1 (p. 33) [7]. Instru-
mentation with longitudinal and transverse strain gauges is applied as in the previous tests of 0° Tension
Test (p. 31).

(5.1)

where EY is the modulus calculated from the initial linear part of the recorded stress-strain rela-
tionship which in principal will have similar form to that shown in Figure 5.3: Typical Recorded
Uniaxial Stress-Strain Relationship (p. 33).

5.1.2. Out-of-Plane Shear

The out-of-plane shear behavior can be characterized by short beam shear experiments. Tests such
as these, using either a three, four or five-point bending set-up, often suffer from the fact that signi-
ficant bending is observed in the sample in addition to the shear deformation.

A test set-up developed at EMI, and used successfully in [2], minimizes bending effects resulting in
maximized shear deformation (see Figure 5.4: Short Beam Shear Test Setup (p. 34)). The shear force
is recorded along with the total displacement using a clip gauge. The shear stress τ and the shear
angle γ can then be derived. For full details of this experiment, see [2].

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 33
Material Characterization Tests

Figure 5.4: Short Beam Shear Test Setup

Figure 5.5: Typical Shear Stress - Shear Angle Relationship (p. 34) shows a typical the shear stress–shear
strain relationship. The out of plane shear stiffness G13 and the matrix yielding stress τmatrix is derived
from the extensometer signals.

Figure 5.5: Typical Shear Stress - Shear Angle Relationship

5.2. Equation of State Properties: Inverse Flyer Plate Tests

Inverse planar plate experiments are used to obtain the uniaxial compression behavior at strain rates
up to 104/s. This method is very well adapted to characterize the shock compression behavior of unknown
materials. Studies have shown the importance of equation of state properties for accurate simulations
of composite materials subjected to hypervelocity impact [1].

A schematic of the test is given in Figure 5.6: Configuration of Inverse Planar Impact Experiments and
a Typical Velocity Trace from the Rear Surface of the Witness Plate (p. 35). A projectile consisting of
cylindrical samples of the composite being tested, covered typically with an aluminum backing, are
accelerated in a gas gun to velocities up to 1000 m/s. The projectile then impacts a stationary witness

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
34 of ANSYS, Inc. and its subsidiaries and affiliates.
 Delamination Properties

plate of well characterized steel. The velocity of the rear surface of the witness plate is recorded using
a high resolution VISAR laser interferometer.

Figure 5.6: Configuration of Inverse Planar Impact Experiments and a Typical Velocity Trace from
the Rear Surface of the Witness Plate

Figure 5.6: Configuration of Inverse Planar Impact Experiments and a Typical Velocity Trace from the
Rear Surface of the Witness Plate (p. 35) shows a typical velocity trace for an inverse flyer plate experi-
ment. After shocking the arriving material sample and the stationary witness plate on the contact surface
to their respective Hugoniot states, the stress waves travel through both plates. The pressure wave inside
the witness plate is converted to a pressure release wave at the free surface and propagates back into
the steel. Due to the impedance mismatch between steel and the sample material, this wave is only
partially transmitted across the impact surface into the sample. It is mainly converted back to a pressure
wave, which provides a further velocity increase when reaching the free surface of the witness plate.
The repeated reflection of the wave at the surfaces results in stepwise increase of the free surface velocity.

The complete impact conditions in the sample can be calculated without any assumption on the mater-
ial’s equation of state. Shock states are determined using the Rankine-Hugoniot equations, and the free
surface velocity can be used to obtain the particle velocity, since the measured free surface velocity ufs
is approximately equal to twice the particle velocity up of the target at the first shock state of Fig-
ure 5.6: Configuration of Inverse Planar Impact Experiments and a Typical Velocity Trace from the Rear
Surface of the Witness Plate (p. 35). Together with the known shock properties of the target, the
Rankine Hugoniot equations can be used to deduce all variables of the shock state. Further details of
this experiment are given in [1][2].

5.3. Delamination Properties

This section describes an overview of experimental techniques used to characterize delamination
thresholds and energies. A number of delamination modes are considered and are illustrated schemat-
ically in Figure 5.7: Delamination Modes. (a) Mode I : Normal Delamination. (b) Mode II : Shear
Delamination. (c) Mode III : Shear Delamination (p. 36).

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 35
Material Characterization Tests

Figure 5.7: Delamination Modes. (a) Mode I : Normal Delamination. (b) Mode II : Shear Delamination.
(c) Mode III : Shear Delamination

5.3.1. Direct Plate Impact - Mode I Delamination Strength

The through thickness strength at high strain rates can be determined by the use of a direct plate
impact test. In this test the projectile and sample are designed in such a way as to ensure that release
waves from the target free surface and from the projectile rear side in the sample (target) material
superimpose resulting in tensile stresses in the sample. If the tensile stress exceeds the loading capacity
of the material, the dynamic through thickness strength under uniaxial strain, or spall strength can
be determined. Mode I delamination occurs and compressive release waves are generated, which
limit the velocity decrease at the free target surface.

Figure 5.8: Experimental Configuration of the Direct Plate Impact Experiment (p. 36) shows the exper-
imental configuration. A cylindrical projectile is accelerated in a gas gun and normally impacts a sta-
tionary sample. A high resolution VISAR laser interferometer records the velocity increase on the target
rear surface.

Figure 5.8: Experimental Configuration of the Direct Plate Impact Experiment

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
36 of ANSYS, Inc. and its subsidiaries and affiliates.
 Delamination Properties

The spall strength can then be calculated from the following equation,
(5.2)

where cp is the soundspeed and Δusp is the difference in velocity between the shocked state and
the limit of spall strength signal.

Figure 5.9: Typical Velocity Trace from a Direct Plate Impact Test

5.3.2. Double Cantilever Beam – Mode I Delamination Energy

The fracture energy for static mode I delamination can be measured using a Double Cantilever Beam,
or DCB test [2]. The test sample has a crack placed into it during manufacture. This predefined crack
is then extended by applying a load to the faces of the sample both above and below the crack. This
can be achieved by gluing hinges onto each surface. These hinges are then clamped into experimental
apparatus which pulls them apart. The test machine used to apply the load to the sample also allows
the sample to move uniaxially along its length by a set of rollers as illustrated in Figure 5.10: Illustration
of the Double Cantilever Beam Test (p. 38).

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 37
Material Characterization Tests

Figure 5.10: Illustration of the Double Cantilever Beam Test

The applied force and the displacement of the cantilevers are used to derive the energy consumed
for crack propagation. The fracture energy GIC is deduced by normalizing over the delaminated area.

A typical force/displacement curve is displayed in Figure 5.11: Force – Displacement Curve for Determ-
ination of the Mode I Fracture Energy Release Rate (p. 38). Before the onset of crack propagation the
slope of the curve is linear. Subsequently a decreasing force level is observed caused by an increasing
lever of the test machine.

Figure 5.11: Force – Displacement Curve for Determination of the Mode I Fracture Energy Release
Rate

For full experimental details, see [2].

5.3.3. Double Notch Shear – Mode II Delamination Strength

The interlaminar shear strength (ILSS), the stress st which Mode II delamination occurs, can be
measured using a Double Notch Shear test (DNS) This type of test should be performed according
to standard proposal ASTM D 3846-79 [10].

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
38 of ANSYS, Inc. and its subsidiaries and affiliates.
 Delamination Properties

Figure 5.12: Configuration of short beam shear test sample

Figure 5.12: Configuration of short beam shear test sample (p. 39) shows the geometry of a typical
DNS specimen. Two notches are machined into the sample, one on either side, up to the same lam-
inate layer. The sample is then axially compressed in a servo-hydraulic testing machine inducing shear
loading on the bond area in between the two notches. Measurements of applied load and displacement
are recorded.

Figure 5.13: Typical Output from DNS Test: Interlaminar Shear Strength

Figure 5.13: Typical Output from DNS Test: Interlaminar Shear Strength (p. 39) shows the typical output
from a DNS experiment. The stress is calculated as the maximum load divided by the cross section
between the two notches of the specimen. The interlaminar shear strength corresponds to the max-
imum stress at failure.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 39
Material Characterization Tests

5.3.4. End Notched Flexure - Mode II Delamination Energy

In Mode II delamination, failure occurs due to excessive shear stress. The End Notched Flexure (ENF)
test determines the elastic energy threshold required for the growth of a predefined delamination
for this mode of failure.

Such a test should be performed to the standard proposal EN 6034 [11]. Figure 5.14: Schematic of the
ENF Configuration Used to Determine the Interlaminar Fracture Energy GIIC (p. 40) shows a schematic
of the experimental set-up and geometry of the sample. The sample, containing a predefined crack
is loaded in a test machine in a three point bending configuration. Displacement of the loading device
and applied load is recorded continuously during the tests.

Figure 5.14: Schematic of the ENF Configuration Used to Determine the Interlaminar Fracture
Energy GIIC

The predefined crack grows further as a result of the Mode II loading: specimen bending and the
resultant shear forces at the crack tip. The total fracture toughness energy GIIC is calculated from the
initial crack length and from the critical load P to start the crack, and, the displacement of the test
machine d at onset of crack extension [11]. All remaining coefficients are sample and test rig dimen-
sions. Further detailed information can be found in [2].

(5.3)

The interlaminar fracture toughness energy is the energy per unit plate width that is necessary to
grow an interlaminar crack. A typical force displacement curve is presented in Figure 5.15: Typical
Load Displacement Curve for Determination of the Mode II Interlaminar Fracture Energy (p. 41) and
shows an interruption of the initial slope in the force displacement curve characterizing the critical
load at the onset of crack propagation in the specimen. Afterwards, the reduced bending stiffness
influences the force displacement curve.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
40 of ANSYS, Inc. and its subsidiaries and affiliates.
 Delamination Properties

Figure 5.15: Typical Load Displacement Curve for Determination of the Mode II Interlaminar
Fracture Energy

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 41
 Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
42 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 6: Derivation of Material Properties
In this chapter we discuss derivation of the material parameters required for the various orthotropic
material modelling options from the experiments outlined in Material Characterization Tests (p. 31).

6.1. Constitutive Properties

The following two sections describe the derivation of the constitutive material properties.
6.1.1. Elastic Properties
6.1.2. Plasticity Parameters

6.1.1. Elastic Properties

The input data required for the orthotropic material model in Autodyn can be entered in two ways.
The user has the option to specify the elastic engineering constants directly by selecting Engineering
for the input type, or the stiffness matrix coefficients may be input by choosing Matrix Coefficients.
This option has been implemented since some stiffness matrix coefficients are obtained directly ex-
perimental tests. For example, depending on the application, C11 may be calculated directly from a
uniaxial strain experiment such as the inverse flyer plate tests of Equation of State Properties: Inverse
Flyer Plate Tests (p. 34).

The stiffness matrix coefficients, in terms of the elastic engineering constants, can be calculated from
the following expressions obtained by inverting the compliance matrix Equation 3.6 (p. 13),

(6.1)

where

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 43
Derivation of Material Properties

Table 6.1: Derivation of Orthotropic Material Elastic Properties (p. 44) outlines the experimental tests
used to calculate values for the engineering elastic constants.

Table 6.1: Derivation of Orthotropic Material Elastic Properties

Property Description
E11 Through thickness Youngs Modulus.

Not measured directly. Can be calculated from

to ensure a consistent stiffness matrix, or
can be estimated from inverse flyer plate tests (see
Equation of State Parameters (p. 47)).
E22 In-plane Youngs Modulus.

Calculated from 0° tension tests.

E33 In-plane Youngs Modulus.

Calculated from 90° tension tests.

ν23 In-plane Poissons ratio.

Calculated from 0° tension tests in which strain

gauges applied in 22 & 33 directions.
ν31 Out-of-plane Poissons ratio.

Calculated from 90° tension tests in which strain

gauges applied in 11 & 33 directions.
ν12 Out-of-plane Poissons ratio.

Often unknown, can be calculated from

if ν12 calculated from 0° tension tests in which strain
gauges applied in 11 & 22 directions and E11 known
or estimated.
G23 In-plane shear modulus.

Average of that calculated from 45° tensile tests

using Equation 6.2 (p. 44).
G12 Out-of-plane shear modulus.

Average value calculated from short beam shear

tests.
G31 Out-of-plane shear modulus.

As G12 above.

The in-plane shear modulus G23 is calculated from the following equation [7] where EY is the modulus
measured in the 45° tension test described in 45° Tension Test (p. 33),

(6.2)

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
44 of ANSYS, Inc. and its subsidiaries and affiliates.
 Constitutive Properties

6.1.2. Plasticity Parameters

As with an isotropic yield criteria that includes work hardening, we can use a uniaxial tension test to
obtain the relationship between and . In practice the master curve can be calculated from
any normal or shear stress-strain test. In this example we choose to use a 0° uniaxial tension test.

In this case the quadratic yield function Equation 3.29 (p. 18) reduces to,
(6.3)

and using the definition of effective stress Equation 3.38 (p. 19), the relationship between effective
stress and σ22 is

(6.4)

In general for uniaxial loading the following relationship results:

(6.5)

and, from Equation 3.33 (p. 19) the associated incremental plastic strain components are
(6.6)

On substitution of Equation 6.6 (p. 45) into the expression for incremental effective plastic strain
Equation 3.40 (p. 20) we find that for uniaxial loadings
(6.7)

For pure shear loading the same arguments result in the following relationship for effective stress
and effective plastic strain increment,
(6.8)

(6.9)

where rr = 44, 55, or 66 depending on shear stress components.

Derivation of the plasticity parameters is best illustrated by an example. In the following discussion
results from tension tests on a 0°/90° woven Kevlar-epoxy composite material are presented [2]. In
this case the stress-strain behavior in the 0°, or 22 direction, is considered to define the master
curve. The plasticity parameter a22 in this case is a free parameter and we are at liberty to set its
value simply to 1.0. Equation 6.5 (p. 45) and Equation 6.6 (p. 45) are now used to transform the
measured values of σ22 and ε22 into effective stresses and effective incremental plastic strains thus
determining the required master curve, as shown in Figure 6.1: Derivation of Master Relationship
from Uniaxial Tension Test (p. 46).

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 45
Derivation of Material Properties

Figure 6.1: Derivation of Master Relationship from Uniaxial Tension Test

Kevlar-epoxy is a 0°/90° woven material and is assumed to be transversely isotropic, therefore a22 =
a33 = 1.0. The in-plane plastic Poisson’s ratio was calculated as 0.26 from 0° tension tests where
strain was measured in both in-plane directions. It was calculated from the average gradient of
measured strains immediately after yielding (see Figure 6.2: Longitudinal Versus Transverse Strain
Measured in 0° Tension Tests. Red Line Indicates Region used to Calculate Value for In-Plane Poissons
Ratio. (Picture Courtesy of EMI [2]) (p. 46)).

Figure 6.2: Longitudinal Versus Transverse Strain Measured in 0° Tension Tests. Red Line Indicates
Region used to Calculate Value for In-Plane Poissons Ratio. (Picture Courtesy of EMI [2])

Therefore from Equation 3.35 (p. 19), a23 can now be calculated as follows:

(6.10)

The out-of-plane PPR was estimated as 0.698 from the 0° tension tests where strain was additionally
measured in the through thickness direction.

Therefore from Equation 3.35 (p. 19), a12 , and hence a13 , can now be calculated as follows:

(6.11)

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
46 of ANSYS, Inc. and its subsidiaries and affiliates.
 Equation of State Parameters

(6.12)

The a44 plasticity coefficient was calibrated through simulation of the uniaxial tension test with
loading at +/- 45° to the fibers. Using the plasticity coefficients derived so far in this section the
quadratic yield function is sufficiently described to allow a44 to be modified until the 45° tension test
results were reproduced.

Best agreement with experiment is achieved with a44 = 4.0 (see Figure 6.3: Results from Simulations
of 45° Tension Tests. (Experimental Result Courtesy of EMI [2]) (p. 47)).

Figure 6.3: Results from Simulations of 45° Tension Tests. (Experimental Result Courtesy of EMI
[2])

In the absence of further experimental data it was assumed the out-of-plane shear plasticity coefficients
are equal to the in-plane properties. Therefore,
(6.13)

6.2. Equation of State Parameters

Figure 6.4: Kevlar/Epoxy IFPT, Influence of Shock Effects. (Experimental Result Courtesy of EMI [1]) (p. 48)
compares the results of simulations of inverse flyer plate tests on a Kevlar-epoxy composite material to
experiment. The blue velocity trace is recorded in a simulation that used an orthotropic model with a
linear equation of state. It is obvious that the initial Hugoniot states are under predicted. However a
much improved correspondence with experiment is achieved when using the orthotropic material
model a polynomial equation of state.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 47
Derivation of Material Properties

Figure 6.4: Kevlar/Epoxy IFPT, Influence of Shock Effects. (Experimental Result Courtesy of EMI
[1])

Based on the results of Figure 6.4: Kevlar/Epoxy IFPT, Influence of Shock Effects. (Experimental Result
Courtesy of EMI [1]) (p. 48) a plot of the derived Hugoniot impact states is given in terms of shock-
versus particle-velocity points in Figure 6.5: Shock Velocity versus Particle Velocity Relationship for Kevlar-
epoxy Inverse Flyer Plate Tests. (Data courtesy of EMI [1]) (p. 48). Further details of this can be found
in [1][2]. Kevlar-epoxy exhibits typical solid behavior with shock velocity approximately increasing linearly
with particle velocity.

Figure 6.5: Shock Velocity versus Particle Velocity Relationship for Kevlar-epoxy Inverse Flyer
Plate Tests. (Data courtesy of EMI [1])

If the application being modelled requires use of a nonlinear equation of state in conjunction with an
orthotropic material model then the gradient and intercept of the line in Figure 6.5: Shock Velocity
versus Particle Velocity Relationship for Kevlar-epoxy Inverse Flyer Plate Tests. (Data courtesy of EMI
[1]) (p. 48) can be directly input as the S1 and C1 inputs, respectively, for a shock equation of state.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
48 of ANSYS, Inc. and its subsidiaries and affiliates.
 Failure and Softening Properties

Alternatively, a polynomial equation of state may be used. In this case the A1 term is calculated as the
effective bulk modulus. The A2 and A3 terms may then be calibrated by simulation to give best agree-
ment with experiment. This is indeed the approach used to obtain the results of Figure 6.4: Kevlar/Epoxy
IFPT, Influence of Shock Effects. (Experimental Result Courtesy of EMI [1]) (p. 48).

Assuming that the 11-direction is defined as being through the thickness of the composite material it
is possible to calculate the C11 stiffness matrix coefficient. Since the flyer plate tests are approximately
uniaxial in strain the following expression applies,
(6.14)

6.3. Failure and Softening Properties

The following sections describe the methods used to determine the failure and softening properties of
materials.
6.3.1. Failure Properties
6.3.2. Softening Properties

6.3.1. Failure Properties

The failure initiation stresses can be derived from experiment as described in Table 6.2: Failure Prop-
erties and the Tests Used to Measured Them (p. 49).

Table 6.2: Failure Properties and the Tests Used to Measured Them

Property Experimental Test

σ11fail Calibrated to experiment to give good agreement
with spall plate tests and DCB.
σ22fail 0° or 90° tension test depending on orientation of
material directions.
σ33fail 0° or 90° tension test depending on orientation of
material directions.
σ31fail Interlaminar shear strength test.
σ23fail 45° tension test and calibration via simulation.
σ12fail Interlaminar shear strength test.

6.3.2. Softening Properties

The fracture energies used to control the rate of crack growth in each cell are obtained from the ex-
periments detailed in Delamination Properties (p. 35) as described below in Table 6.3: Fracture Energies
and Tests From Which They May be Derived (p. 49).

Table 6.3: Fracture Energies and Tests From Which They May be Derived

Property Experimental Test

Mode I delamination 11 - Gf11 DCB test.
Gf22 Calibrate simulation to 0° tension
test.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 49
Derivation of Material Properties

Property Experimental Test

Gf33 Calibrate simulation to 0° tension
test.
Mode II delamination 13 – Gf13 ENF test.
Gf12 (same as 13 direction) ENF test.
Gf23 Assumed to be same as other
shear directions.

6.4. Calculating Laminate Properties from Uni-Directional Data

Laminates made from composite materials often consist of a number of repeating sublaminates or indi-
vidual lamina. Typically material data may be available for a uni-directional lamina and it is therefore
desirable to be able to be able to calculate the effective material properties for the laminate from the
constituent lamina properties.

This is also useful since within Autodyn each layer of the laminate is not explicitly modelled, rather
continuum elements representing equivalent homogeneous anisotropic solids are used to represent
thick laminates consisting of a number of repeating lamina.

The approach of Sun and Li [4] has successfully been used to calculate laminate properties for use in
Autodyn and is now described.

Consider a laminate consisting of N orthotropic fiber composite lamina of arbitrary fiber orientations.
In the following description the x and y co-ordinates are in the plane of the composite and z is through
the thickness. The effective macro-stress and macro-strains are defined to give,

(6.15)

(6.16)

where and are the stresses and strains in the kth lamina and, if t k is the lamina thickness and
h the total thickness
(6.17)

The effective elastic properties for the laminate are given by,

(6.18)

The x-y plane for a lamina is a plane of symmetry and is also a symmetry plane for the effective solid.
Therefore the effective stiffness matrix reduces to the following form,

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
50 of ANSYS, Inc. and its subsidiaries and affiliates.
 Calculating Laminate Properties from Uni-Directional Data

(6.19)

and effective compliance matrix is calculated as,

(6.20)

where the effective elastic properties are then,

(6.21)

After lengthy algebraic manipulations the following expressions are recovered for the effective stiffness
matrix coefficients [4],

(6.22)

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 51
Derivation of Material Properties

where,

The principal material directions of all the lamina in the laminate will not necessarily be aligned with
the global axes. Therefore, before the above summations can be performed it will be necessary to
transform the stiffness matrix for each lamina as follows,
(6.23)

where the transformation matrix T is defined as,

(6.24)

and the angle θ is the angle of the material directions with regards to the global axes.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
52 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 7: Example Applications
The composite material models available in Autodyn and described throughout this report have been
validated through application to a wide variety of loading conditions. In this section of the report a
selection of these applications are presented.

7.1. Hypervelocity Impacts

7.1.1. Advanced Material Model for Hypervelocity Impact Simulations

This summary describes the work carried out by the Ernst-Mach-Institute (EMI), Century Dynamics Ltd
(CDL) and TNO Prins Maurits Laboratorium (TNO) under ESA research contract No. 12400/97/NL/PA(SC)
[1]. The aim of this project, termed AMMHIS, was to develop and define advanced material models
and data for Nextel and Kevlar/Epoxy under hypervelocity impact conditions. The AMMHIS project
was particularly concerned with the materials for the spacecraft shielding configuration used to protect
the European Columbus module of the ISS. Due to the experimentally observed behavior of Nextel
and Kevlar-epoxy, a new orthotropic hydrocode model was developed. This was necessary because
pressure inside these materials depends on deviatoric strain components as well as volumetric strain.
Nonlinear effects, such as shock effects, can be incorporated through the volumetric straining in the
material. The developed model includes orthotropic material stiffness and a nonlinear equation of
state and is described in Equations of State (p. 14).

The hypervelocity impact of aluminum projectiles up to 15mm in diameter and normal velocities in
the range of 3 km/s to 15 km/s were considered. The quality of the model has been demonstrated
by comparison of simulations with tests to characterize the involved materials as well as with impact
tests on the complete shielding configuration.

Figure 7.1: Inverse Flyer Plate Tests, Experimental and AMMHIS Model Results

The developed material model is able to predict the main aspects of the shielding material response.
Calculated shielding damage in the first bumper and in the Nextel and Kevlar-epoxy layers correlates
well with experimental results. In terms of the back wall damage; for the 3km/s impact the simulation
predicts a small hole that was not observed in the experiment, for the 6km/s impact the predicted

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 53
Example Applications

back wall damage is consistent with the experimental observations (Figure 7.2: Alenia/EMI Test A8611
– Material Status During Impact on Reference Shielding, 15mm Diameter Projectile, 6.5km/s (p. 54)
and Figure 7.3: Alenia/EMI Test A8611- Key Features of Material Response During Impact on Reference
Shielding, 15mm Diameter Projectile, 6.5km/s (p. 54)). The simulation sensitivity studies carried out
at both these velocities suggest that these are marginal cases in terms of back wall penetration.

Figure 7.2: Alenia/EMI Test A8611 – Material Status During Impact on Reference Shielding,
15mm Diameter Projectile, 6.5km/s

Figure 7.3: Alenia/EMI Test A8611- Key Features of Material Response During Impact on Reference
Shielding, 15mm Diameter Projectile, 6.5km/s

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
54 of ANSYS, Inc. and its subsidiaries and affiliates.
 Hypervelocity Impacts

7.1.2. Advanced Material Damage Models

This summary describes the work carried out by the Ernst-Mach-Institute (EMI) and Century Dynamics
Ltd (CDL) under ESA research contract No. 12400/97/NL/PA(SC), CCN No. 2 [2]. The aim of this project,
termed ADAMMO, was to develop and further improve the composite material characterisation and
modelling techniques for modelling hypervelocity impact conditions based on the study of Advanced
Material Model for Hypervelocity Impact Simulations (p. 53). The focus of the study was to allow
better prediction of the levels of damage and residual strength in composite materials after impact.
A methodology was implemented to incorporate anisotropic hardening observed in some highly
nonlinear composite materials such as Kevlar-epoxy, whilst retaining the ability to couple the con-
stitutive response to a nonlinear equation of state. Simulations of static tensile tests agreed well with
experiment (see Figure 7.4: Results from Simulation of 0° Tension Test (p. 55)).

Figure 7.4: Results from Simulation of 0° Tension Test

An orthotropic damage model, outlined in Orthotropic Damage Model (p. 28), was implemented into
Autodyn during the ADAMMO project. In short the softening behavior observed in some composite
materials is modelled using a crack softening approach where damage is accumulated as a function
of crack strain and the ability of a laminate to carry tensile loads is reduced.

A simulation was performed with similar loading to that experienced by the material in a short beam
bending test.

Figure 7.5: Simulation of Short Beam Shear Test

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 55
Example Applications

The resulting shear stress-deflection curve obtained from the simulation is shown in Figure 7.5: Simu-
lation of Short Beam Shear Test (p. 55). The response is linear elastic until the material begins to fail
in shear after which a much reduced stiffness is initially observed. As the shear strain in the material
continues to increase, the response hardens as the fibers re-orientate and start to pick up the shear
load. This behavior is the same as that observed in the short beam bending characterisation test and
demonstrates the versatility of the orthotropic model.

A number of the plate impact damage tests were simulated. In this test a symmetric assembly of
Kevlar-epoxy and Al plates are used to prevent delamination caused by superposition of release
waves. A schematic diagram of the test geometry is shown in Figure 7.6: Schematic Diagram of Exper-
imental Setup (p. 56).

Figure 7.6: Schematic Diagram of Experimental Setup

Figure 7.7: Simulation Results of 276 m/s Impact Velocity

(a) Orthotropic damage (b) Orthotropic damage (c) Brittle damage model
model model (AMMHIS)

Through thickness Through thickness plastic

damage strain

The orthotropic damage model predicts some damage for a 276m/s velocity impact, due to hardening
of the material (Figure 7.7: Simulation Results of 276 m/s Impact Velocity (p. 56)). This is close to the
experimental result that shows limited amounts of delamination. The brittle damage model shows
extensive delamination and even some bulk failure of the material which is clearly an over prediction
of the levels of damage.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
56 of ANSYS, Inc. and its subsidiaries and affiliates.
 Hypervelocity Impacts

In conjunction with the hardening model described above, these additional orthotropic material
models have been extensively validated by comparison with hypervelocity impact experiments. A
series of debris cloud damage experiments were performed at EMI designed to generate impacted
targets with damage gradients that are fully damaged in the central impact region and little damage
at the extent of the targets.

Figure 7.8: Test 4355 – Through Thickness Damage During Impact (p. 57) shows the simulation of
test 4355 - a 8.2mm Al sphere impacting an 2.0mm Al bumper at 4.68km/s. The resulting debris cloud
then impacts a 5.7mm thick Kevlar-epoxy target.

Figure 7.9: Final Damage of Test 4355 (p. 58) shows the final state of the simulation and the same
configuration simulated without hardening included and using the brittle damage model (AMMHIS).
The ADAMMO simulation produces much more compact level of damage with significant amounts
of intact material remaining adjacent to the main impact zone as observed in experiment.

Figure 7.8: Test 4355 – Through Thickness Damage During Impact

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 57
Example Applications

Figure 7.9: Final Damage of Test 4355

a) 4355 b) T04355 - ADAMMO c) BM4355 - AMMHIS

7.2. Ballistic Impact Examples

Composite and textile armor systems are increasingly being utilized as impact protection materials in
weight critical environments. Applications range from protective shielding for space vehicles against
hypervelocity impacts, to personal protective equipment of the soldier against ballistic threat. The use
of composite and textile armor systems can result in a reduction in weight while maintaining impact
performance, or increased impact performance for a given weight. The limitations of performing con-
trolled impact experiments, especially at higher velocities, on composite and textile materials means
that there is a need to study these events using numerical simulations.

Numerical investigations of ballistic fragment impacts on Polyethylene based fiber composite (Dyneema)
[12] and Kevlar/Epoxy [13] body armor have been performed at Century Dynamics and are briefly de-
scribed below and compared with experimental work.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
58 of ANSYS, Inc. and its subsidiaries and affiliates.
 Ballistic Impact Examples

7.2.1. Ballistic Fragment Impacts on Aramid Composite Plates

An orthotropic material model coupled with a polynomial equation of state has been used to conduct
numerical investigations of FSP (Fragment Simulating Projectile) impacts on Kevlar fiber Reinforced
Panel (KFRP) systems used for personal protective clothing armor [13]. Results of these investigations
are shown below and compared with instrumented ballistic impact experiments carried out at TNO
Prins Maurits Laboratory. The experiments included measurement of target back surface velocity
during the impact using VISAR techniques.

In this work, the coupled anisotropic material model was used to simulate the case of a 1.1g steel
fragment impacting an aramid composite plate at 483m/s. The Lagrange processor of Autodyn-2D
was used to represent both fragment and composite target.

The aramid composite is made up of 19 layers of woven Kevlar 29 bonded in Epoxy. Basic quasi-
static elastic properties for the composite were known although no dynamic material properties that
would allow the derivation of nonlinear shock terms were available. In the initial simulation work re-
ported in [13], the material properties for Kevlar-129/Epoxy derived in [1] were used along with the
polynomial equation of state. The 4340 steel fragment was represented using a Johnson-Cook strength
model and data.

Figure 7.10: Fragment Impact on KFRP at 483m/s, Simulation and Experimental Results

a) Simulated Damage Development b) Back surface Velocity

7.2.2. Polyethylene Fiber-Based Armor

The performance of the orthotropic model coupled with a nonlinear equation of state has been further
evaluated for the ballistic performance of a new Dyneema (UD-HB25) composite, which was being
assessed by DCTA and TNO through an experimental testing programme [12]. Simulations of inverse
flyer plate tests were conducted and a remarkable correspondence between test and simulation was
obtained using the new composite model. The derived model was then successfully used to perform

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 59
Example Applications

simulations of fragment impact tests and make an assessment of the V50 impact velocity. The model
also reproduced the deformation and delamination extent of the target plate, as illustrated in Fig-
ure 7.11: Autodyn Simulations of 1.1g FSP Impacting 3.2mm Dyneema UDHB25 (Magenta Regions
Indicate Delamination) (p. 60).

Figure 7.11: Autodyn Simulations of 1.1g FSP Impacting 3.2mm Dyneema UDHB25 (Magenta
Regions Indicate Delamination)

7.3. Bird Strike Example using Composite Shell Elements

A layered composite shell has been implemented in Autodyn-3D. Composite shell elements allow efficient
modelling of thin composite materials that are subjected to structural (rather than shock) type loading.
All of the existing material models that are compatible with the standard shell element can be applied
to individual layers of new composite shell elements. In addition the orthotropic equation of state can
be used with the layered composite shell element. The composite shell always assumes that the 11 and
22 directions are in the plane of the shell and the 33 direction lies through the thickness.

An example composite shell application is shown in Figure 7.12: Initial Configuration of Composite Tail
Section (p. 61) to Figure 7.14: Material Status Following Bird Strike (p. 61). Results from an Autodyn
simulation of a bird strike onto the leading edge of a composite tail plane are shown. Figure 7.12: Initial
Configuration of Composite Tail Section (p. 61) and Figure 7.13: Initial Configuration of Composite Tail
Section: Impact Zone (p. 61) show the initial configuration of the bird and tail section. The bird has
been modelled with SPH particles and the tail plane with a combination of composite and standard
shell elements. The composite shells have been used to represent the layered GFRP and CFRP skin of
the tail section. Figure 7.14: Material Status Following Bird Strike (p. 61) shows material status plots for
the tail section following the impact. Results for the outer and inner layers of the skin are shown.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
60 of ANSYS, Inc. and its subsidiaries and affiliates.
 Bird Strike Example using Composite Shell Elements

Figure 7.12: Initial Configuration of Composite Tail Section

Figure 7.13: Initial Configuration of Composite Tail Section: Impact Zone

Figure 7.14: Material Status Following Bird Strike

(a) Outer sublayer (b) Inner sublayer

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 61
 Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
62 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 8: Recommendations
Autodyn contains many options for the modelling of composite materials. These options are described
in detail in this report. However, in using these complex models the problems are often two fold;
knowing which options/models to use and how/where to obtain the required material properties.

In Orthotropic Constitutive Models (p. 11), the various models are described in detail. It is important
to take a pragmatic view of them rather than simply selecting the most complex. For example, if the
application being modelled is only subjected to a relatively low speed impact then using an orthotropic
material model with a linear equation of state is sufficient. In this case relatively simple experiments
are required to obtain the directional strength properties or manufactures material data may be sufficient.

If the application is at ballistic velocities or higher then shock effects are most likely to be important
and using a nonlinear equation of state is recommended. However, this requires additional experiments
as described in Equation of State Properties: Inverse Flyer Plate Tests (p. 34), in addition to the direc-
tional strength properties.

Whether or not the hardening option is required depends only upon the material being modeled. Some
materials, such as carbon-epoxy composites have been observed to exhibit a very linear behavior [14].
Other composite materials, like Kevlar-epoxy, are highly nonlinear and significant hardening is observed
[2]. In this case it is imperative to have the stress-strain data obtained from the experiments of Direc-
tional Strength Properties (p. 31) so that the plasticity parameters can be calculated as outlined in
Plasticity Parameters (p. 45). Since the tension tests are usually performed to ultimate failure, this ex-
perimental data will also indicate whether or not softening is significant.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 63
 Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
64 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 9: References
1. S. Hiermaier, W. Riedel, C.J. Hayhurst, R.A. Clegg, C.M. Wentzel, Advanced Material Models for Hypervelocity
Impact Simulations, EMI-Report No. E43/99, ESA CR(P) 4305, 1999.

2. W. Riedel, W. Harwick, D.M. White, R.A. Clegg, Advanced Material Damage Models for Numerical Simulation
Codes, EMI-Report No. I 75/03, ESA CR(P) 4397, October 2003.

3. J.K. Chen, [Link], [Link],“A Quadratic Yield Function for Fiber-Reinforced Composites“, Journal of
Composite Materials, vol. 31, No. 8, 1997.

4. C.T. Sun and S. Li,“Three-Dimensional Effective Elastic Constants for Thick Laminates”, Journal of Composite
Materials, vol. 22, 1988.

5. C.E. Anderson, P.A. Cox, et. al.“A Constitutive Formulation for Anisotropic Materials Suitable for Wave
Propagation Computer programme-II”, Comp. Mech., vol. 15, p201-223, 1994.

6. J.P. Hou, N. Petrinic, C. Ruiz,“Prediction of Impact Damage in Composite Plates”, Comp. In Sc. and Tech.,
60, 273-281, 2000.

7. R.M. Jones,“Mechanics of Composite Materials”, McGraw-Hill, 1975.

8. M.A. Meyers,“Dynamic Behavior of Materials”, John Wiley & Sons, 1994 – ISBN 0-471-58262-X.

9. W.C. Kim, C.K.H. Dharan, Analysis of five-point bending for determination of the interlaminar shear strength
of unidirectional composite materials, Composite Structures, Vol. 30, 1995, 241-251.

10. ASTM standard D3846-79 (1985), ASTM Standards and Literature References for Composite Materials, 2nd
Ed., American Society for Testing and Materials, Philadelphia, PA, 1990.

11. Luft- und Raumfahrt, Kohlenstoffaserverstärkte Kunstoffe, Prüfverfahren, Bestimmung der interlaminaren
Energiefreisetzungsrate, Mode II,G/IIC, DIN EN 6034, April 1996.

12. C.J. Hayhurst, J.G. Leahy, M. Deutekom, M. Jacobs, P. Kelly,“Development of Material Models for Numerical
Simulation of Ballistic Impact onto Polyethylene Fibrous Armour”, presented at PASS, Colchester, UK, 2000.

13. R.A. Clegg, C.J. Hayhurst, J.G. Leahy, M. Deutekom,.“Application of a Coupled Anisotropic Material Model
to High Velocity Impact of Composite Textile Armour”, presented at the 18th Int. Symposium on Ballistics,
San Antonio, USA, 1999

14. D.M. White, E.A. Taylor, R.A. Clegg,“Numerical Simulation and Experimental Characterisation of Direct Hy-
pervelocity Impact on a Spacecraft Hybrid Carbon Fiber/Kevlar Composite Structure”, Int. J. Impact Eng.,
vol. 29, 2003, pp779-790.

Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 65
 Release 2020 R1 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
66 of ANSYS, Inc. and its subsidiaries and affiliates.