Finite element modeling of static indentation damage in laminated composites


Material Information

Finite element modeling of static indentation damage in laminated composites
Physical Description:
xi, 134 leaves : ill. ; 29 cm.
Jung, Han Soo, 1952-
Publication Date:


bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1991.
Includes bibliographical references (leaves 108-111).
Statement of Responsibility:
by Han Soo Jung.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 001693388
notis - AJA5467
oclc - 25223151
System ID:

Full Text








This thesis is dedicated to my father-in-law, Dong Seung

Kim, who could not wait to see this thesis. He passed away in

August 1989.

Also this thesis is dedicated to my younger brother-in-

law, strong and courageous seafarer, Captain Ja Sup Kim, who

must have perished with his crew sometime during September

1990 on his last voyage from Chile.


First of all, acknowledgement must go to the professors

of my supervisory committee. Without the care and guidance of

Dr. Bhavani V. Sankar, chairman of the committee, this work

could not have been possible. Dr. Lawrence E. Malvern, who

also was a member of the supervisory committee for my master

program, encouraged me to continue my graduate study to the

doctoral program. Dr. C. T. Sun introduced me to composite

materials, about which I had no experience before. Dr. Ibrahim

Ebcioglu taught me about plate theories. His encouragement

always helped me when I was dispirited. Dr. John J. Mecholsky

has widened my vision about composite materials. His many

suggestions greatly helped this research. I also heartily

appreciate the financial assistance I received during the

doctoral program.

My special appreciation extends to Dr. C. A. Ross and Dr.

J. E. Milton at the Graduate Center of the University of

Florida at the Eglin Air Force Base in Fort Walton Beach,

Florida. Both professors were members of my supervisory

committee for the master's degree. Without their help, I could

not have begun my graduate education. So many fond memories

are there with them.


I would like to express my appreciation for the financial

support of the United States Air Force for my master's

program. I appreciate Dr. David Ebeoglu, Mr. Robert J. Arnold

and Dr. David J. Butler of the Armament Division, Eglin AFB,

Florida. Dr. Ebeouglu was always there whenever I need help.

Mr. Arnold was very kind to me and to my family. The advice

and cheerfulness of Dr. Butler helped me all the time.

I sincerely appreciate the Agency for Defense Development

where I devoted fourteen years of my young life. I extend

special thanks to Dr. J. H. Hong, Dr. Y. K. Yoon, Dr. Y. H.

Kang, Dr. K. S. Chon, Dr. S. H. Do and many other bosses and

colleagues: I wish I could name them all.

Finally, I am very grateful to my wife Ki Young and two

children, En Jae and En Kyum, who must hold the traveling

mileage records in their schools, for their endurance and

support. And special appreciations must go to my parents,

brothers and sisters for their prayers and encouragement.



ACKNOWLEDGEMENTS ....................................... iii

LIST OF TABLES ............. .............................. vii

LIST OF FIGURES ....................................... viii

ABSTRACT ................................................. x


1 INTRODUCTION .......... ........................ .. 1

1.1 Introduction ...................................... 1
1.2 Experimental Background ........................ 5
1.3 Objective and Scope ............................. 7

CONTACT PROCESS ................................ 15

2.1 Introduction .... ...................... ........ 15
2.2 Axisymmetric Formulation .................... 17
2.3 Constitutive Relations for an Axisymmetric
Laminated Plate ............................... 19
2.4 Algorithm for the Contact Program .............. 25
2.5 Axisymmetric Finite Element ................... 28


3.1 Introduction ............................... 44
3.2 Evaluation of the Program ................... 45
3.3 Comparison with Test Results ................ 48
3.4 Comparison to Contact Laws .................. 49
3.4.1 Indentation Law .. ...................... 50
3.4.2 Contact Pressure Distribution .......... 54
3.4.3 Contact Radius ........................... 54
3.4.4 Center Deflection ....................... 55
3.5 Friction Effect ............................ 57

4 DAMAGE ANALYSIS ............................... 78

4.1 Failure Criteria ............................... 78

4.2 Implementing Delamination into Finite
Element Program ................................ 88
4.2.1 Introduction ............................ 88
4.2.2 Numerical Experiment ................... 89
4.3 Implementing Fiber Failure and Matrix Failure
into Finite Element Program ................... 90
4.4 Incremental Finite Element Method for
Damage Analysis ................................ 91

5 CONCLUSIONS .................................. 106

5.1 Conclusion ................................... 106
5.2 Comments for the Future Research .............. 107

REFERENCES ............................................ 108

APPENDIX ................................................. 112

BIOGRAPHICAL SKETCH .................................... 134



2.1 Material properties of AS4/3501-6 ........... 41

3.1 Material properties of the plates
used in. the examples ....................... 60

3.2 Material properties of the
numerical applications ...................... 67

4.1 Ultimate strength of AS4/3501-6 lamina ...... 87



1.1 Impact test scheme ............................ 10

1.2 Test configuration ............................ 11

1.3 Contact force vs. center deflection curve ...... 12

1.4 C-scan picture of impacted composite plate ..... 13

1.5 A photo-micrograph of the impacted plate ....... 14

2.1 The continuity of the stress through the
thickness of the composite plate ............... 38

2.2 The continuity of the strain through the
thickness of the composite plate ............... 39

2.3 Stacking of plies to one element ............... 40

2.4 Indentation of a plate ........................ 42

2.5 Axisymmetric element ........................... 43

3.1 Plate dimensions in numerical examples ......... 59

3.2 Contact force vs. contact radius relationship .. 61

3.3 Interlaminar shear stresses in the
graphite/epoxy laminates ....................... 62

3.4 Interlaminar shear stresses
in the hybrid laminates ........................ 63

3.5 Contact force vs. center deflection curve ...... 64

3.6 Contact force vs. center deflection curve,
elastic range ................................. 65

3.7 Configurations for numerical examples .......... 66

3.8 Contact force vs. indentation curve,
isotropic material ............................. 68


3.9 Contact force vs. indentation curve,
composite plate (4 mm thickness) ............... 69

3.10 Effect of plate radius on contact force vs.
center deflection curve ....................... 70

3.11 Effect of plate radius on contact force vs.
indentation curve ............................. 71

3.12 Effect of indenter radius on contact
force vs. indentation curve .................... 72

3.13 Contact pressure distribution .................. 73

3.14 Contact force vs. contact radius ............... 74

3.15 Contact force vs. contact radius, comparison ... 75

3.16 Contact force vs. center deflection............. 76

3.17 Friction effect ................................ 77

4.1 Partial failure patterns ...................... 86

4.2 Delamination simulation ....................... 97

4.3 Numerical application of
delamination effect ........................... 98

4.4 Delamination effect on contact force vs. center
deflection curve ............................... 99

4.5 Geometry of the contact surface ............... 100

4.6 Contact force vs. center deflection curve ..... 101

4.7 Contact force vs. center deflection curve ..... 102

4.8 Damage pattern of the composite plate ......... 103

4.9 Damage pattern of the composite plate ......... 104

4.10 Damage pattern of the composite plate ......... 105

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy





Chairman : Bhavani V. Sankar
Major Department: Aerospace Engineering,
Mechanics and Engineering Science

The objective of this research is to study progressive

failure of circular laminated brittle matrix composite plates

due to static indentation loading using finite elements. The

plates are assumed to be simply supported along the edge, and

are indented by a rigid indenter with a hemispherical nose. A

quasi-isotropic homogenization process is used to make the

problem axisymmetric. Four noded, isoparametric,

quadrilateral, axisymmetric solid elements are used to model

the plate. A new algorithm is used to analyze the progressive

contact between the indenter and the plate. Hashin's failure

criteria are modified so that they can be applied to

axisymmetric, transversely isotropic laminates. Stiffness

matrices of failed elements are modified suitably. New slip

surface is introduced when delamination occurs. Very stiff

springs are attached to nodes on delamination surfaces to

prevent interpenetration of the nodes. The results from

incremental damage analysis agree well with experimental

results. This program can be applied to analyze impact

problems involving large masses impacting at very low




1.1 Introduction

Currently composite materials are used in primary

structures of aircraft, particularly aircraft wings and

frames, because significant weight savings can be achieved due

to their high strength to weight ratio. Also, composite

materials can easily be tailored to obtain desired properties

in different directions. For example, more 00 layers can be

used for longitudinal loads, 900 layers for transverse loads

and 450 layers for shear loads. Also, there are other

advantages such as resistance to damage by fatigue loading,

immunity to corrosion, ease of processing, and low coefficient

of thermal expansion. Structural forms that are otherwise

inconvenient or impossible to manufacture, e.g., more complex

airfoils or helicopter blades, can be easily fabricated.

Hence, the utilization of composite materials in aerospace and

transportation industries is continuously increasing [1,2].

However, it is well known that brittle matrix composites

such as graphite epoxy are very susceptible to low velocity

impact, which may be caused by dropping tools, runway debris,

hailstones, bird strikes and ground service vehicles hitting


aircraft structures. High and medium levels of impact energy

cause surface damage that is relatively easily detected, and

therefore, repairs may be undertaken. However, low velocity

impacts which are at quite low energy levels can seriously

damage laminated composite structures internally so that

frequently it is very hard to detect damage through visual

inspections. The damage can be detected only by using

nondestructive inspection (NDI) techniques. For example, low

velocity impact to an aluminum plate may cause denting damage

and the structure does not lose very much strength. However,

inside of an impacted composite plate there may be fiber

breakage, matrix cracks, and delaminations, while the impacted

surface may not show any damage. Commonly it is termed as

"Barely Visible Impact Damage" or BVID [2]. After impact,

composite materials lose considerable strength because of

internal damage and delaminations. [2-7]

It becomes very important to understand the response of

composites to the impact of rigid objects and to predict

damage in order to make composite materials more reliable. The

study of the effect of BVID on static and fatigue strengths of

graphite/epoxy panels shows that impact degradation is of more

concern in compression loading than in tension because of

delaminations which occur during impact. Importance of the

impact study is emphasized by the fact that upper surfaces of

the aircraft are more likely to receive impacts from tool


drops than lower surfaces, and upper surfaces operate in a

compression-dominated load regime [2].

Because of the extremely complicated nature of impact on

the composite materials, empirical methods have been widely

employed. Generally, empirical methods have several distinct

limitations. First, it is very expensive and also time

consuming to conduct testing and to process the data. Second,

it is practically difficult to test all the combinations of

material properties, boundary conditions of the composite

structures, and impactor velocities and weight variations.

This problem can be critical, if the composite material is in

a developing stage, and only limited quantities are available.

Third, it is very difficult to comprehend all the stress and

strain distributions, and damage progression inside composites

by empirical methods. And also there are difficulties in

matching results from specimens to the real structures. [6,7]

Therefore, numerous analytical methods have been

developed to overcome the shortcomings of empirical methods.

However, so far, because of the extremely complicated nature

of impact problem, there are no complete three-dimensional

models available which can simulate the entire impact process

and predict damage progression. There is an urgent need to

develop a three dimensional progressive finite element program

which calculates contact forces, contact displacements,

stresses and strains inside the composite plate. During the

contact process, the finite element program must be able to


check all the failure criteria for fiber breakage, matrix

cracking and delamination. If failure occurs, then the program

could modify the stiffness matrix to account for failures .

Fortunately, low velocity impact due to large masses can

be treated as a static indentation problem, because the impact

duration is much longer than the time required by the

propagating waves to travel from the impact site to the

supports or free edges, and the test results show good

agreement with this assumption. Strictly speaking, the

terminology of low velocity impact is not clearly defined.

Rather low impact velocity is generally defined as the impact

velocity which is low enough to justify a static analysis of

the response of the structure, especially in the vicinity of

impact [3,4,8,9]. This velocity depends upon the mass, shape

and the material property of the impactor and also upon the

impacted structure boundary conditions and its natural

frequencies. Generally it ranges 1-3 meter/second with

impactor mass less than 10 kg [10].

The advantage of static analysis is that the time

variable is eliminated. Hence, we can closely look into

progressive damage mechanisms with a small amount of

computational resources. If time effects are considered, it

will be a formidable task to obtain all the meaningful stress

and strain distributions inside the plate and to map the

damage pattern and its progression because of computational

time and storage limitations of the computer system. If


coarser elements are employed, it will be very hard to obtain

accurate data because meaningful data must be obtained in the

small confined area under the indenter.

The experimentally determined contact laws account for

local damage in an empirical manner, and hence cannot be used

to study the local damage. A fully three dimensional

formulation is required to study the progressive damage.

In the present study a finite element program performs an

incremental analysis of the contact problem, and assesses the

damage due to indentation.

1.2 Experimental Background

In this section some results from an experimental study

[11] are presented to bring out the importance and necessity

of the finite element modeling of the damage process due to

static indentation.

Test specimens were made of Hercules carbon prepreg tape

AS4/3501-6, which is an amine-cured epoxy resin reinforced

with unidirectional carbon fibers. All specimens were cut from

a 30.48 cm (12 inch) by 30.48 cm, [00, +450, -450, 90]8s

composite plate lay-up. In order to obtain the specimens with

the quality suggested by the manufacturer, the curing

instructions was carefully followed as in Ref. [12]. An

autoclave (Baron Model BAC-24) was used to fabricate all the

specimens of composite laminates used in this study. Detailed


procedure of fabrication can be found in Ref. [11]. Material

property of the 00 lamina can be found in Table 2.1.

Figure 1.1 shows the pendulum impact test facility in

which all impact tests were conducted. The pendulum includes

the support strings connecting the tup to the laboratory

ceiling, and the tup, which is a combination of weights, a

load cell, and an indenter with one inch diameter

hemispherical head. The total mass of the tup is 13.84 kg.

Static indentation tests were done in a standard MTS material

testing machine.

The composite laminate was placed on a circular steel

ring of 5.08 cm (2 inch) diameter to simulate a circular plate

with simply supported boundary condition. The test

configuration is shown in Figure 1.2.

The test results of contact force vs. center deflection

of the circular composite plate for the static indentation and

impact tests are shown in Figure 1.3. This graph suggests that

static indentation tests can simulate low velocity impact

tests. Also it suggests that indenting process is highly

nonlinear due to internal damages. A sudden drop of the

contact force is a unique feature of indentation of composite

plates. Figure 1.4 shows the picture of the C-scan of the

impacted plate. The circular damage pattern justifies the

axisymmetric formulation used in the present study. Figure 1.5

shows the picture of the photo-micrograph of the polished


cross section of the indented plate. This shows matrix cracks

due to shear failure and delaminations.

However, these data are not sufficient to understand the

entire impact or indentation process. They do not explain how

damage initiates and propagates. Also, it is very hard to

quantify this damage.

1.3 Objective and Scope

The major objective of this research is to develop a

finite element program which models the static response of

laminated composite plates indented by a rigid spherical

indenter and to predict progression of damages inside the

plate. The results are applicable to low velocity impact of

the composite plates.

A circular plate and a rigid hemispherical indenter were

chosen in modeling because they correspond to many practical

impact situations. Quadrilateral, isoparametric, axisymmetric

finite elements were employed in this research. By employing

axisymmetric elements, the three dimensional problem can be

treated as two dimensional by eliminating dependency on 6.

In Chapter 2, the finite element formulation and the

incremental displacement method are introduced. The input to

the finite element program is the displacements necessary to

close the gap between the indenter and candidate contact nodes

of the plate. The laminated plate consists of many plies that

are orthotropic and have different angles of fiber


orientation. As axisymmetric finite elements were employed in

this research, there must be special considerations of the

material properties and failure criteria to account for this

discrepancy. The in-plane strains, exx, e and yx are

continuous at the interface between two adjacent plies whereas

the stresses, a a and are not. On the other hand, zz,

cxz and yz are continuous through thickness and are better

represented in a homogeneous plate model than ezz, Yxz and yyz,

which can be discontinuous through the thickness. A set of

transversely isotropic elastic constants is derived for a

group of plies that can be considered as a quasi-isotropic


In Chapter 3, the accuracy and efficiency of the present

method were evaluated by comparing the results with some early

work [10]. The contact force vs. center deflection curve of a

AS4/3501-6 graphite/epoxy plate was compared with experimental

results in the elastic range. Then, some widely used contact

laws were compared with the present finite element results.

The effect of friction was studied by assuming stick friction

conditions between the indenter and the top surface of the


In Chapter 4, the failure criteria are introduced. Among

many interactive and independent failure modes, three-

dimensional failure criteria introduced by Hashin [13,14] were

employed in the present studies. These criteria have the same

quadratic form as Tsai-Hill's criterion, but they can

distinguish different failure modes such as fiber breakage,

matrix cracking and delamination. Since we use an axisymmetric

formulation, the damage pattern is also going to be

axisymmetric, which does not strictly represent the damage

pattern in a composite laminate. However, experimental impact

studies [15] have indicated that the actual damage pattern is

like a spiralling stair case, and the projection of these

delaminations can be considered as circular and hence

axisymmetric. This was also confirmed by studies in the

laboratory [11]. The extent of delamination, fiber breakage

and matrix failure were mapped during the contact process.

After failures are detected, the total stiffness matrix of the

plate must be modified properly depending upon the failure

modes. Matrix and fiber failure can be implemented by reducing

or eliminating element stiffness matrices. Implementing

delamination into the finite element program requires creation

of new free surfaces, consequently new nodes. Rigid springs

are inserted across the newly created free surfaces to prevent

interpenetration of the nodes. Delamination can be assumed as

Mode II (Shear Mode) crack propagation. Incremental damage

analysis results are compared with experimental results.

A summary of conclusions and some suggestions for future

research are listed in Chapter 5.

Figure 1.1 Impact test Scheme


E 2
e o

....... ...... I .......... ... .

R = 12.7 mm
7 4 mm

Figure 1.2 Test configuration

0 2 4 6
Center Deflection (mm)

Figure 1.3 Contact Force vs. center deflection Curve

Figure 1.4 C-Scan picture of impacted composite plate

Figure 1.5 A photo-micrograph of the impacted plate



2.1 Introduction

The finite element method has been established as a

powerful and versatile tool of analysis for contact problems

involving isotropic bodies. In the present study, a method for

solving the problem of indentation of composite laminates is


The contact problem of quadratic surfaces of elastic

bodies was first solved analytically by Hertz in 1881. Since

then, numerous contact problems have been solved by

formulating the problem as an integral equation. The solution

is either obtained by solving the equation numerically or, for

simple configurations, in closed form. The Hertzian type

problem has three general assumptions. First, the contacting

solids are homogeneous, isotropic and linear elastic. Second,

the contact areas are essentially flat and relatively small

compared to the various radii of curvature of the undeformed

bodies in the vicinity of the contact interfaces. Third, the

contacting solids are perfectly smooth and only normal

pressures are considered [6,16].


Hertzian type solutions are not suitable for indentation

of composite laminates which are inhomogeneous, anisotropic

and also involve significant damages at low load levels.

Further, most laminated composites in use can not be

adequately represented as a half-space as in Hertzian contact

problem. Because of these complexities, three dimensional

analysis, especially three dimensional finite element methods,

are required in the present study.

In the last two decades, many researchers thought that

contact problems could be considered as a special case of

constrained minimization of either total or complementary

potential energy. The minimization is formulated as a

mathematical programming problem, and the solutions are

obtained by either using incremental linear programming,

quadratic programming or conjugate (modified) gradient

methods. The minimization problem can also be formulated as a

variational inequality using penalty methods or Lagrange

multipliers [17-21]. However, these methods have certain

disadvantages such as creation of flexibility matrix, the

awkwardness involved in introducing additional unknowns, such

as slack variables and Lagrange parameters. Nothing can be

better than direct formulation of finite element method as a

simple design tool, which calculates the detailed stress field

in various types of laminates for various boundary conditions

accounting for the damages inside the plate. The goal is to

use simple and conventional elements so that the method can be


implemented in commercial finite element programs. The novel

incremental contact algorithm developed in this study is

explained in Section 2.4.

2.2 Axisymmetric Formulation

The benefit of employing an axisymmetric element is that

the three dimensional problem can be treated as two

dimensional. Especially in the case of rigid ball indentation,

where the contact surface is circular, axisymmetric

formulation is a very efficient method to treat the indenting

process in a computer program, because upper surface nodes of

the plate come into contact successively and the contact area

is readily determined. If we use three dimensional solid

elements, it will be a formidable task to determine the

contact nodes for the case of a spherical indenter, although

the employment of three dimensional solid elements has a

distinct advantage in that we can use orthotropic material

properties and also the failure criteria for orthotropic

laminates directly without any modifications. Complexity of

the contact algorithm demands employment of axisymmetric

elements, and consequently there must be some considerations

to convert the material properties of a laminated plate, which

are basically orthotropic, into equivalent axisymmetric

properties. The circular damage distribution pattern of the C-

scan pictures of the damaged plates (Figure 1.4) justifies the

employment of axisymmetric elements to solve the problem.


Stress-strain relationships of the axisymmetric element are

discussed extensively in Section 2.3.

In the present study u(r,z) and w(r,z) represent the

radial and transverse displacements, respectively. The

axisymmetric stress and strain components are

{ o } = { ra zz 'rz ee } (2.1)

(e} ={ ez ezz Yrz eee } (2.2)


err- r

zz W (2.3)

au aw
Y +
v = 8z + Br
a" z ar

considered as elastic until failure, material non-linearities

As graphite/epoxy composite is very brittle and can be

considered as elastic until failure, material non-linearities

are ignored [2]. As will be seen later, the deflections

considered are smaller than the plate thickness, and hence

geometric non-linearities are also not considered.


2.3 Constitutive Relations for an Axisymmetric Laminated Plate

A composite laminate can be considered as having several

sublaminates. If each sublaminate is quasi-isotropic, then

each sublaminate can be idealized as a homogeneous

transversely isotropic material. Then the composite laminate

can be considered as a laminate made up of several

transversely isotropic layers. Hence the axisymmetric

formulation can be employed if the external loading and

boundary conditions are also axisymmetric.

The stress-strain relationship in the Cartesian

coordinates are

{e} =[S]{o (2.4)

where ( o } and ( e ) are,

{ e } = { e eyy Y~ ezz Yzx Y~ } T (2.5)

{ a } = { Oxx ayy, xy Ozz Tzx rzy }T (2.6)

For the ply which has fiber orientation angle 8 from the

x coordinate, lamina constitutive relations can be obtained


{e} = [ T][S][ T {o }e


where [ S ] is

1 V21 0 -v3 0 0

V12 1 0 320 0
E1 E2 E3
0 0 0 0 0

0 o o
3 230 1 0 0
0 0 0 0 0
o o 0 0 0

and [ To ], [ Te ] are given by,

cos20 sin28 2cos6sinO 0 0 0
sin28 cos2 -2cosisinO 0 0 0
-cos6sinO cosOsinO cos28-sin2e 0 0 0
0 0 0 1 0 0
0 0 0 0 cos6 sin6
0 0 0 0 -sine cos6

cos2) sin28 cosfsin6 0 0 0
sin28 cos28 -cosssine 0 0 0
-2cos6sin6 2cos6sinO cos28-sin28 0 0 0
0 0 0 1 0 0
0 0 0 0 cosO sine
0 0 0 0 -sin6 cosO

For the fiber orientation angle 6,

{ e = [ S]e {o }e

[ S ]e = ] [ s ]e[ S]T ]T

Averaging [ S ]e for the number of plies from this form

is not suitable because in-plane stresses such as axx, ao and

Txy are not continuous through thickness of the element while
exx, e and Yy are continuous. And, also, out-of-plane

stresses ozz, zx and Tzy are continuous through thickness while

out-of-plane strains ezz, yzx and yzy are not continuous. These

are well illustrated in Figure 2.1 and Figure 2.2, which Wu

and Springer have obtained [23]. These figures show the

continuity of the stresses and strains through the thickness

of the impacted plate. The best way to average the modulus of

the element can be done after denoting [ S ], as,

P q
[S], = (2.8)
r s

where p, q, r and s are 3 x 3 matrices. Then the stress-strain

relationship can be written as

(xx exx
ayy yy

x p q (2.9)

oz [ r s ezz
tzx zx
tzy zy

The stress-strain relationship can be rearranged such that all

the components of stress and strain that are continuous

through the thickness are on one side and the remaining

components are on the other side as shown in Figure 2.9.

ax e
a e

ezz rp- s-rp-q zz

Szy zy

Let us denote [ D ]e

[-D ], = (2.11)
Srp- s-rp-qj

Then an average [ D ] matrix for a group of N plies can

be obtained as,

[D]av= [D]e,


Now the elastic matrix [ C ], can be recovered from

[ D ]av. The [ C ]av can be expressed in terms of P, Q, R and

S as,

P-QS-1R QS-1 (2.13)
[ e S- R S-1

The stress-strain relations take the form

{ 0 }av = [ C ] { e }a

For widely used cross ply or quasi-isotropic laminates,

such as [ 0 90 ]s or [ 00 ,+ 450 900 ]i laminates, [ C ]av

becomes transversely isotropic. For axisymmetric problems,


are = 16e = re = e3 = 0



"Crr Ce 0 Crz

C]e Ce Cee 0 C (2.15)
0 0 C o
cz C, 0 Cz

Further, the present homogenization scheme yielded a

symmetric stiffness matrix, i.e., Cij = Cji, and also

c'r = ce

Crz = Cez (2.16)

Sometimes, five engineering elastic constants, such as

Er, Ez, Grz, Vre, Vrz can describe the material properties of a

transversely isotropic laminate plate [10]. For this case, the

stress-strain relationship can be obtained as follows.

1 v 1
1 Vre 0
0 -Z

Vre 1 0 Vzz
E, E E
[C]=r E E r (2.17)

V V 1
0z rz 0
Er E, E,
r r


Table 2.1 shows the material properties of the AS4/3501-6

graphite/epoxy lamina. Also, the calculated transversely

isotropic element material property matrix [ C ],, for a plate,

which consists of four plies in [ 00 450 900 ] orientation

is given in Table 2.1. This configuration is chosen because of

its wide application in aircraft structures. For example, the

main wing box for a straight wing aircraft, 00 plies carry the

spanwise direct stresses induced by wing bending like a

cantilever beam bending, 450 plies carry shear stresses

caused by wing twisting, and finally 900 plies carry the

chordwise direct stresses due to fore and aft bending of the

aircraft wing [2].

2.4 Algorithm for the Contact Problem

Unlike conventional finite element problems, contact

problems are mixed boundary value problems where nonzero

displacements are specified at parts of the boundary and zero

tractions are specified in the rest. Further, the

displacements in the contact surface have to satisfy certain

conditions specified by the geometry of the indenter. Sankar

[10] used a contact algorithm in conjunction with a finite

difference method. In the present study, the algorithm has

been modified to suit the finite element method. Unlike the

finite difference method, the finite element method has the

advantage of modeling complex geometries and multi-layer



For the frictionless contact problem, the vertical

degrees of freedom of the top surface nodes of the plate

become contact candidates. In Equation (2.18), the

displacements of these contact candidate degrees of freedom

are notated as du and other displacements as dr. The dr degrees

of freedom can be condensed to obtain a smaller stiffness

matrix as shown below.

[K. Kur du Fu (2.18)
.ru Krr dr 0

{ d } = [ Kr ]-1 [ Kru ] { du } (2.19)

After condensation, the equilibrium equations take the form

[ Kc ] { d = { Fu } (2.20)


[ KC ] = [ K,] [ K ] [ K ]-[ Kru (2.21)

Since we are dealing with an axisymmetric problem, the

contact area is always circular. Let us assume that the

contact has progressed to Node k-1 (see Figure 2.4). The


procedure for determining the incremental load required to

bring the kth Node into contact is as follows. Equation (2.20)

can be partitioned as

K.= Kct w. Fe (2.22)
Kfe Kff, wf 0 J

where wc are the z displacements of nodes within the contact

region including k-1 and wf are the remaining degrees of

freedom for which the forces are zero. Let us use gi to denote

the gap between the indenter and nodes outside the contact

region, i.e., i 2 k. The first step is to determine the

incremental displacement of the indenter to close the gap gk.

We first apply a unit displacement to the indenter, i.e.,

Wi = 1, i = 1,2, ,k-l (2.23)

The equation (2.22) is solved for Fc and wf. Now the new gap

g'i becomes

g'i = gi 1 + w., i t k (2.24)

Now we can calculate the displacement 6 of the indenter

to make the gap g'k identically equal to zero.

6= gkk
gk g k (2.25)

1 Wk )

After determining 8, the nodal force FC can be computed. By

integrating all the nodal forces, we obtain the incremental

contact force, AP. The procedure is repeated for the next

contact node k+l.

When frictional effects are included, radial forces at

nodes in the contact region can not be considered as zero. The

procedure has to be slightly modified, as is explained in

Section 3.5.

The damage alters the stiffness matrix. We assume that

the indentation process is under displacement control, and

hence the w-displacements of the contact nodes are kept the

same when passing from one iteration to the next. The changes

in [ K ] cause a change in the contact load (usually

reduction) for the same indentation displacement. Detailed

damage analyses are presented in Chapter 4.

2.5 Axisymmetric Finite Element

In this section, the derivation of the stiffness matrix

of a quadrilateral, isoparametric, axisymmetric element is

given for the sake of completeness.


In isoparametric elements, the element coordinates and

the element displacements are expressed with the same shape

functions associated with nodal coordinates and displacements

respectively as shown below.

r=C N= (Ej') **


n (2.26)
u = N (,1) ui
w = Ni (in) w_

where r and z are element coordinates, & and n are the natural

coordinates and u and w are element displacements as shown in

Figure 2.5. Subscript i signifies node number. Shape function

or interpolation function Ni are expressed in the natural

coordinate system, which has variable & and t that each vary

from -1 to +1 as shown in Figure 2.5. The characteristic of

the shape function Ni is that its value is unity at node i and

is zero at all other nodes. For the quadrilateral element, the

Ni are as follows.

IN, = ( 1 + E ) ( 1 + r )

N2 = ( 1 ) ( +

V = (1 +- ) (1- )

To evaluate the stiffness matrix of the element, we need

to calculate the strain-displacement transformation matrix for

the small strain case as,

err = -ar


Su aw

ee =

As the interpolation functions, N,, are defined in the

natural coordinate system by Equation (2.27), we need the

relationship between r and z derivatives and t and q

derivatives. By using chain rules, the following expressions

can be obtained.

au zar z1 au au
I a aj tW I= -1 at (2.28)
au ar az au au
z 9 z] owI ar}

= as, asF an = as (2.29)
aw ar az aw aw
BT =W RW -1 a1 (2.29)
aw r az aw aw
Fz _F _T1 anrl

In the above equations, J is the Jacobian operator
relating the natural coordinate derivatives to the local
coordinate derivatives. The inverse of Jacobian matrix must
exist in order for the isoparametric element to be used. This
means that there is a one-to-one correspondence between the
natural and the local coordinates of the elements. In actual
computations, the sign of the Jacobian determinant is
monitored at the Gaussian integration points to ensure the
existence of the inverse Jacobian. Generally, a zero or
negative value of the Jacobian determinant indicates input
data error or an overly distorted element, and the
computations are terminated automatically for the purpose of
node renumbering [22].
Equations (2.28) and (2.29) can be calculated as follows.

ar ( 1 1 1
aE 4 4 4 4
r _1( 1 1 1
ar 1(l+)r1- (I-F) r2- (l-) r3-- (l+-) r,
4 4 4 4

1 (1+ ) z- (1+1) Z2- (1-n) z3+-! (1-1) Z4
z 1 (+)-1 1 +1 -

47 4 4 4 4

u= _1 (+()u- 1 1 u+
au +4 1 (1, -u1- (1+ ) u2 --1 (1-1)U3

(1-n) u4

au 1 (1+4) u1 + (1-&) u2- (l-) u3- (+ ) U4
4 4 4 4


aw = (1+ ) w- 1 (1+11 2 (1- )3
a 4 4 4

(1-11) w4

1 1 1
aw= (1+) w+ (1-) (1- W3- (1 ) w4
-l 4 4 4 4

Eventually, the unknown variables of the system will be

the nodal displacements { U }

{ U} = { UI w 2 W2 u3 w3 u4 w4 }T



f 9U. ,

u 4 1+ o 0 i- o -(-+) o -((1+- o0 0
au 7 + 0 o- 0 (1-d) 0 (1+n) 0 (-U

r P_ 1j 01 + 1 (1+l) 0 -(1-) 0 l-rT
aw 4 0 1+t 0 1- 0 -(1-t) 0 -(1+)

Strain can be expressed with nodal variables { U ), and
also stress can be obtained by multiplying the strain vector
by the elasticity matrix [ C ]e

{ e = [B] { U}
{ o } = [ ]e [B]e { U}

In the axisymmetric case, unlike other isoparametric
elements, the strain-displacement matrix [ B ]e has terms
containing r. Thus, for two elements of identical shapes and
material properties but at different r locations, the [ K ]e
will be different. Also, a problem arises when r is zero
because some terms have r in the denominator. A proper
numerical scheme, such as Gaussian quadrature, can avoid this
problem [24].
In the axisymmetric case, { e ) becomes

{ e } = {err ezz Yrz eee}T (2.36)

and [ B ]e becomes 4 by 8 matrix as follows.

1+8 0 1++
0 1i+ 1+1 0
(0 1 & (1-&) (1+11)
-(1+ ) o i- -+
1 0 1-( -(i+n) 0 (2.37)
4 _(i_ ) 0 -( -S) -1 ) ( + )
o -(1-S) -(1-a) 0
1-n 0 -(1+0)
0 -(1+) 1-r 0

The element stiffness matrix can be formulated by using

the virtual work principle. The principle requires that for a

body under equilibrium,

f 8eTodV = f 8aTfBdV+ f a'fsdS + 8aTfN (2.38)
vs v, S,


a ; actual stress

be ; virtual strain that corresponds to the imposed

virtual displacement

8a ; virtual displacement

fa ; actual body force

fs ; actual surface force
f ; actual nodal force

for the axisymmetric element, the volume integral of any

function A becomes

f A dV= f A rdrdzdO (2.39)

Volume integration can be written in the natural

coordinate system. Notice that r, which is a global

coordinate, appears in the integration, which is a unique

feature of the axisymmetric element.

f A d = 2nff A r det [J] dd4 (2.40)

If there are no body forces, or initial stresses,

Equation (2.38) becomes

6aT [ B ]e[ C ][ B ]e dV { U} = aT fX (2.41)

by applying

{ 8e } = [B ]e { U}

{ ao = [ C]e { = [ C ] [ C] [B ] { U} (2.42)
{6a } = [ N] {BU}

The stiffness matrix is obtained as,

[K] = [B ][ C][B ]e dV (2.43)

It is very easy and convenient to use Gauss-Legendre

numerical quadrature integration formula to evaluate Equation

(2.43). The Gauss-Legendre integration points Ei and ij and

corresponding weighing factor aij are found in Ref. [24]. Care

must be taken to substitute r in Equations (2.37) and (2.43).

From Equation (2.26), r can be obtained as

1 1
r = (1+Ei ) (1+Il)r. +- (1-+) (1+,).r2
*+ (i-S) (1- .)i r3+- (1+i) (-j ) r4
4 4

After substituting Ei and mj into E and at the previous

equations, [ K ]e can be obtained numerically as


[K ]e = [ B ]i[ Ce[ B aijdet[ J ]i,r
i,j = 1

The assembly of element stiffness matrices, applying

boundary conditions and the solution of resulting linear

equations follow the standard procedures that can be found in

Ref. [24,25].


.. -20
.0.04 0 0.04



Figure 2.1 The continuity of the stresses through the
thickness of the composite plate, Ref. [23]

0 0.04







.0.04 0 0.04


O.OOS i------

-O.OOS 1-'- -; ---- .- ,






-0.04 0 0.04

x, COORDINATE ( in )

Figure 2.2 The continuity of the strains through the
thickness of the composite plate, Ref. [23]



One Element Consists of Four Plies

0 Degree

-45 Degree

45 Degree

90 Degree

Laminated Plate

Figure 2.3 Stacking of plies to one element

Table 2.1 Material Properties of AS4/3501-6

(unit; Gpa)

E, 144.8 E2=E3 9.7

G12=G13 6.0 G23 3.6

V12= 13 0.3 23 0.34

* Courtesy of Hercules Advanced Materials and

Systems Company

Axisymmetric Element Property Matrix,

for [ 00, 450, 900 ]

[C ],

(unit ; Gpa)

63.593 4.178

4.178 11.105

0 20.198

0 4.178

0 4.500

4.178 0 63.593



a; Radius of the plate
c; Contact radius
R; radius of the indenter
Gk; Gap between indenter


k- c k- k IGk N

Figure 2.4 Indentation of a plate

.h.( ,r, u

One Element

2 =1 1 3 4

1 (= 1

2 1
3 =-1

Local Coordinate z Global Coordinate

Axisymmetric Element

Figure 2.5



3.1 Introduction

In this chapter, the accuracy of the finite element

program is evaluated by comparing with the results of Sankar

[10], who employed a finite difference method to solve the

equilibrium equations. The results obtained using only 384

finite elements (24 x 16) in the present study compared very

well with that of Ref.[10]. In Section 3.2, available methods

of increasing toughness and damage tolerance are briefly

discussed. In Section 3.3, the finite element program was

applied to the experimental test configuration shown in Figure

1.2 to obtain contact force vs. center deflection. Without

modification of the total stiffness of the plate due to

internal damages, finite element results seem to deviate from

the experimental results. These deviations can be attributed

to various damages as will be discussed in Chapter 4. However,

in the elastic range, where damage is minimal, finite element

results match the test curves well (Figure 3.6). In Section

3.4, various empirical equations including the Hertzian

indentation law are briefly introduced and compared with

finite element results in the elastic range. In Section 3.4.1,


the finite element program is used to study the indentation of

solid cylinders and plates to compare with Hertzian

indentation law. Effects of plate radius and indenter radius

on indentation and deflection of the plate are studied. In

Section 3.4.2, Hertzian contact pressure distribution, which

is known to be elliptic, was compared with finite element

results. In Section 3.4.3, the relation between contact force

and contact radius was studied. It is known that in a half-

space contact force varies as the cube of the contact radius.

Finite element results matched very well with the available

empirical equations. In Section 3.4.4, the relation between

contact force and center deflection was studied. In Section

3.5, the effect of friction between indenter and plate was


3.2 Evaluation of the Program

Hybridization, such as mixing the fiber types, thus

utilizing the advantageous properties of one material to

overcome deficiencies of another, is one method to deal with

impact damages. The higher strain to failure of glass fiber or

Aramid fibers (Kevlar) led to substantial increases in

threshold energies for impact damage, with the energy level

required to produce delamination being increased. The

Aramid/graphite/epoxy hybridization is promising because

Aramid material is relatively inexpensive and it increases

impact resistance while graphite fibers provide compressive


strength, which Aramid fibers lack. Tests showed that doubling

of residual strength after impact was possible using

Aramid/graphite hybridization compared to an all graphite

specimen [2].

It is well known that tough resins can significantly

reduce the damage caused by impact and substantially improve

the residual strength following impact. This is because the

onset of delamination and the delamination growth in a

composite depends on its interlaminar fracture toughness,

which in turn depends on the toughness of the matrix material.

Widely used epoxy resins are brittle, but crystalline

thermoplastics such as PEEK (Poly-ether-ether-ketone) can be

used in place of epoxy to increase toughness. Mode I fracture

toughness of the graphite/epoxy ranges 80-200 J/m2 while the

fracture toughness of graphite/PEEK is about 2500 J/m2 [2]

Interleaving is a process in which thin layers of

thermoplastics are placed in between layers of thermosetting

prepreg materials to form a tough interlayer. It is useful if

the failure occurs through delamination. It is not effective

if the transverse shear is high and, for this case, an

improvement of shear modulus of the resin is more important


Through thickness reinforcements, such as a three

dimensional braiding process or through-the-thickness

stitching, can improve impact resistance. However, in-plane


strength will be reduced because of wavy fibers and irregular

fiber orientations.

Sankar [10] performed the analysis of indentation of

hybrid and interleaved plate problems by using a finite

difference method starting from the axisymmetric equilibrium


C11 ( ,rr + r-1U, r r-2 ) + C44 U,z
+ ( C13 + C44 ) W, = 0

( C13 + C44 ) ( u,Irz + r-1u,z ) + C33W, z

+ C44 ( w,I + r-1W, r ) = 0




{rr C11
(aee = C12
ozz C13

C13 rr
C13 ee
C33 ezz

trz = C44 Yrz


The results obtained by the present finite element

program coincide very well with results in Ref. [10]. Figure

3.2 shows contact force vs. contact radius relationship for



the three different laminated plates. Material properties of

the composite laminates are listed in Table 3.1.

Figure 3.3 shows interlaminar shear stress in the

homogeneous graphite/epoxy laminate. Figure 3.4 shows

interlaminar shear stress in the hybrid laminate. Employing

the finite element program is very economical because only 384

(24 x 16) elements were used, and computing time was about 30

minutes in a Vax main frame computer. Although this analysis

does not account for the internal damages of the plate, it

provides information about the possible onset of delamination

and its subsequent propagation. As delamination failure occurs

where shear stress is high, shear stress distributions such as

Figures 3.3 and 3.4 suggest locations where delamination

failures possibly may initiate. By obtaining other strain and

stress distributions inside plate, possible damage pattern,

such as matrix crack and fiber breakage, and its extent can be

predicted conservatively and less expensively before detailed

damage analysis. Strengthening by inserting tough materials or

hybrid laminates to minimize delaminations can be studied by

using the finite element program.

Detail damage analyses are presented in Chapter 4.

3.3 Comparison with Test Results

The finite element program was used to analyze some of

the specimens used in the indentation tests discussed in

Section 1.2. The contact force vs. center deflection curve was


obtained without considering damages inside the plate and was

compared with the test results in Figure 3.5. Even in the

presumably elastic range, finite element results make the

plate appear stiffer. However, from Hashin's failure criteria,

which will be discussed in Section 4.1, the damages begin to

accumulate at about 140 N contact force or 0.33 mm center

deflection. Hence the actual elastic range is very small.

Comparing the results only in this elastic range, results of

finite element and experimental data showed fairly good

agreement. A total of 512 elements (32 x 16) were employed.

Fine mesh was used near the center of the plate, and coarser

elements were employed outside the contact region. Actually,

finite element results form an upper bound for the test

results. It can be easily explained that these early

deviations of the stiffness of the test results are due to the

inherent imperfections in the test specimens such as voids,

debonds and nonuniform fiber orientations.

3.4 Comparison of Contact Behavior

In this section, various results concerning indentation,

contact pressure distribution, contact radius and center

deflection obtained using the present finite element method

are compared with available results.

3.4.1 Hertzian contact law

Because of the extremely difficult nature of analytical

solutions of contact problem, which involves permanent

deformations and damage such as matrix cracks and

delamination, an experimental approach has been taken to

determine the contact law of composite materials. Naturally,

loading and unloading curves (contact force-indentation

relation) are different because damages occur even at low load

levels. This is because composite materials, unlike monolithic

materials, have very low resistance to the damage due to loads

transverse to the fiber direction.

For the loading curve of contact between an elastic

sphere and elastic half space, the Hertzian contact law [6] is

F =K a2 (3.5)

where Kc is

4( 1 vs )2 ( 1 v )2
Kc I[ E + E1- (3.6)

Subscript s denotes the indenter and t denotes the target. In

Equation (3.6), RS, v and Es are the radius, the Poisson's

ratio and the Young's modulus of the indenter, respectively.

When a metal indenter like a steel ball indents a polymer


composite, the indenter can be assumed rigid because Young's

modulus of steel is much higher than Ez, the modulus of

elasticity of the composite in the transverse direction. For

loading, Sun and Yang [6] suggested the contact force-

indentation relationship in the form

= a 2 (3.7)

A more complex form of contact law obtained by Conway

[26] for an elastic spherical impactor striking a transversely

isotropic composite medium has only 3 percent difference from

Equation (3.7) in the example problem to be presented in this

section. In Equation (3.7), C44 of [ C ] matrix can be

replaced by Ez without significant error (4 percent difference

in this example problem).

The unloading formula suggested by Sun and Yang [6] was

F= c (a-ao)25 (3.8)


a0 = a [1-( acr)] am > ac
am (3.9)

a = 0 am acr


In Equation (3.9), am is the indentation corresponding to the

maximum contact force before unloading. It is assumed that the

critical indentation, ac, depends only on the material


Contact force vs. indentation graphs were obtained by the

finite element program for the clamped solid cylinder of 20 mm

radius and 20 mm thickness of a hypothetical isotropic

material, and also for AS4/3501-6 composite laminate plate of

25.4 mm radius. The material properties are listed in Table

3.1. A total of 512 elements (32 x 16) were employed. Inside

the contact region, fine mesh is used, while coarser mesh size

is assigned outside the contact region. Each computer run took

75 minutes on the average in the Vax main frame computer. In

the case of indentation of the thick cylinder, the center

displacement at the top becomes indentation directly. For the

plate, indentation is defined as displacement difference

between top and bottom surface of the plate at the center,

a = w(0) w(h) (3.10)

The results for contact force-indentation are presented

in Figures 3.8 and 3.9. For the solid cylinder of isotropic

material, the finite element results and Hertzian indentation

law of Equation (3.6) agree very well as seen in Figure 3.8,

although errors were expected because we have not solved the

half space contact problem, and the indenter can not be


simulated exactly as a sphere in a finite element program. For

the plate of isotropic material, contact forces predicted by

the finite element model are larger than forces calculated by

Hertzian contact law when the indentation is large (Figure

3.8). Sankar [27] noticed this difference in the contact

coefficients of Equation (3.5) between contacts of a sphere

with a half space and with a plate. He suggested that

differences of the definitions of indentation between half

space and plate may be the cause of the problem. The finite

element contact forces of the composite plate (4 mm thickness)

are shown in Figure 3.9.

To check the influence of the radius of the plate, 25.4

mm, 38.1 mm, and 50.8 mm radius plates are compared in Figures

3.10 and 3.11 (the indenter radius is 12.7 mm). Although the

contact force vs. center deflection curve shows a big

difference due to the bending effect of the plate (Figure

3.10), contact force vs. indentation curve does not vary so

much (Figure 3.11). At small indentation, contact force seems

to be independent of thickness and radius of the target.

To verify that contact force is proportional to the

square root of indenter radius, two different indenter radii,

6.35 mm and 12.7 mm, were used to indent a 4 mm composite

plate (Figure 3.12). The contact force curve obtained using

the finite element program for the indenter radius 6.35 mm was

denoted as F1 and the curve for 12.7 mm indenter as F2. It may


be seen that for the small indentation, the contact force is

proportional to the square root of indenter radius.

3.4.2 Contact pressure distribution of the plate

The contact pressure distribution for half-space

indentation is as follows [27].

S=3 P (r/C)2 (3.11)
2 3t C2

where P is contact force, and C is contact radius. Finite

element results were compared with these equations in Figure

3.13. Obtaining contact stress distribution by experimental

test is not easy and no data is available for the tests done

as in Section 1.2. Contact pressure distribution calculated

by the finite element method in the 4 mm plate (Figure 3.13)

is elliptic. Contact pressure of the finite element method was

computed from the average vertical stress oz of the elements

in contact with the indenter. Although calculated contact

forces were larger than predicted by the half-space formula,

Equation (3.11), the contact pressure distribution is still


3.4.3 Contact radius of the plate

In a half-space the contact force is proportional to the

cube of the contact radius C. By inserting an empirically


obtained relationship between indentation and contact radius

for the indenter radius as in Equation (3.12) into Hertzian

contact law, Equation (3.7), we can obtain Equation (3.13)


a (3.12)

= z C3 (3.13)

Finite element results for a 4 mm plate are compared with the

above equation in Figure 3.14. In Figure 3.15, contact force

vs. contact radius curves of plates of different thicknesses

are compared. As plate thickness becomes larger, more contact

force is needed to maintain the same radius of contact, as

expected. However, for the small contact radius, contact force

can be assumed to be independent of thickness of the plate as

suggested in Equation (3.13).

3.4.4 Center deflection of the plate

The center deflection 5 of the impacted surface of the

plate can be expressed as the sum of indentation and of

displacements due to bending and shear. In the case of small

deflection, membrane effects can be neglected. Indentation was


obtained from Equation (3.7). Plate bending stiffness, Kp, can

be expressed as [29],

16x (l+v)D


for the axisymmetric plate with thickness h and radius a. In

Equation (3.14), D is


D 12 2)
12 (1-v2)

where, E can be replaced by C1 of the matrix [ C ]av and v,

Poisson ratio, can be replaced by v12. Shear deflection can be

obtained from

aw P
r = G -
zz ar 2trh

Plate shear stiffness, Ks, can be obtained as

. 2x Gh
[ Ln a Ln C]

and C can be obtained from Equation (3.12).

Then, the center deflection 6 is expressed as



S1.5 + P + P (3.18)
KK, K,

Finite element results of a 4 mm plate are compared with

Equation (3.18) in Figure 3.16. It is shown that Equation

(3.18) can predict good results for the elastic response of

the plate. This result suggests that center deflection is a

more reliable parameter than indentation to obtain the contact


3.5 Friction Effect

The effect of friction is studied by assuming that the

coefficient of friction is very large. There are only

compressive forces between the upper surface of the plate and

the indenter. Stick condition requires that there be no

relative motion during a load increment between indenter and

plate at the contact nodes. This condition is satisfied by

constraining the radial displacements of the nodes after they

come into contact with the indenter [30-31]. This can be

easily implemented in the finite element program by modifying

the contact algorithm which was developed in Section 2.4. Now,

displacements in Equation (2.22) in Section 2.4 must include

radial displacements of the upper nodes. In Equation (2.22),

wC is replaced by

{ wk } = { u1, w,, u1 ,2, .... Ui-, W-1} T



Treating the vertical displacements in the same way as in

Chapter 2.4, we now put zero for the radial displacements in

Equation (3.19) during contact.

Radial force is very large compared to vertical force

because of the high modulus of radial and circumferential

stiffness due to the presence of fibers. Hence, matrix failure

may occur virtually immediately after contact starts.

Calculated contact force vs. center deflection for both

frictional and frictionless contact are shown in Figure 3.17.

The effect of friction is significant at higher contact

forces. However, the radial stresses developed may cause

failure even at smaller loads, and the result shown in Figure

3.17 may not be applicable in practical situations.






Rs; Radius of indenter, 10 mm

h; Thickness of plate, 2 mm

C ; Contact radius

a ; Radius of plate, 50 mm

Figure 3.1 Plate dimensions in numerical examples

Table 3.1 Material properties of the plates used in the
examples [10],
(unit, GPa)

Material Er Ez Grz vrz V r

1 50 7 3.8 0.285 0.265

2 20 8 3.6 0.270 0.255

3 3 3 1.1 0.350 0.350

1. Homogeneous laminate plate

material 1, thickness 2 mm

2. Hybrid laminate plate

material 1, thickness 0.25 mm

material 2, thickness 0.25 mm

material 1, thickness 1.5 mm

3. Interleaf laminate plate

material 1, thickness 0.25 mm

material 3, thickness 0.25 mm

material 1, thickness 1.5 mm

3E-01 5E-01 7E-01 9E-01

Contact Radius (mm)

Figure 3.2 Contact force vs. contact radius relationship













150 -
S 140

S 110
L 100 c/h 0.5
L 90
0 80
(0 70 -
To c/h 0,376
60 -
c/h 0.25
40 -
30 -
c/h 0.125
20 -

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Z Coordinate (mm)

Figure 3.3 Interlaminar shear stresses in the graphite/epoxy

170 -
150 -
S, c/h0.5
0 140 -
130 -
P 120 -
S 110
L c/h0.375
S 100 -
L -0
) 70 -
60 -
40 -
30 -
20 -
10 -
0 -
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

Z Coordinate (mm)

Figure 3.4 Interlaminar shear stresses in the hybrid

0 0.2 0.4 0.6 0.8 1
Center Deflection (mm)

Figure 3.5 Contact force vs. center deflection curve

S i

/ r


*.L //






* F.E.M.
..... Impact Test(IMPCB)
....... Static Test(STC1O)
- Static Test(STC 13)
-- Statcl Test(STC 14)

S .08


Center Deflection (mm)

Figure 3.6 Contact force vs. center deflection curve,
elastic range


I /

, -



220 -

200 -

180 -


140 -

120 -

100 -



Rs: Radius of Indenter

^ a ----.\

Solid Cylinder
Rs = 12.7 mm
a = 20 mm


Figure 3.7 Configurations for numerical examples


Table 3.2 Material properties of the Numerical applications

Isotropic Material

E = 2.0 GPa

v = 0.25

G= 0.8

S E(1-v)
S (1+v) (1-2v)

[ C] =





0 1-2v
2 (1-v)

v v
1-v 1-v

0.8 0
2.4 0
0 0.8
0.8 0


Composite Material

Er = 56.427 GPa

Ez = 10.688 GPa

Grz = 4.5 GPa

rz = 0.263

v = 0.3

[C]av= 0


4.178 0 20.198
11.105 0 4.178
0 4.5 0
4.178 0 63.593






140 -

130 0 F.E.M.(4 mm Plate)

S120 0 F.E.M.(Solid Cylinder)
c- B
S 110 -- Eq. 3.6

Z 100 -

90 -
o 80
S 70 -

50 -

40 -

30 -

20 -

0 0.02 0.04
Indentation (mm)

Figure 3.8. Contact force vs. indentation curve, isotropic




130 N F.E.M.

120 Eq. 3.7

O 110

| 100 -

o U

480 -

270 x
,6 80 -

4 0 -E W

10 NW

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016

Indentation (mm)

Figure 3.9. Contact force vs. indentation curve, composite
plate (4 mm thickness)



350 -


z /
0 250 -'
S 200
0 o

0 150 -

--a = 25.4 mm

/ a = 50.8 mm
so Zr

^ ~--- i i i i -- i i i -
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

Center Deflection (mm)

Figure 3.10. Effect of plate radius on contact force vs.
center deflection curve


350 -

S-- a = 25.4 mm / /
o 300 -
---... a = 38.1 mm

5 -- a = 50.8 mm /
260 -
O /.
o //
U. /,"

200 -
0 /SOI



0 ----- I I I ti I
0 0.004 0.008 0.012 0.016 0.02 0.024

Indentation (mm)

Figure 3.11. Effect of plate radius on contact force vs.
indentation curve


350 -

S- F1 (R = 6.35 mm)
3 ........ F2(R= 12.7 mm)
300 /
z F2-F1

S 250
0 /

,( 200 /
r /"



0 0.004 0.008 0.012 0.016 0.02 0.024

Indentation (mm)

Figure 3.12. Effect of indenter radius on contact force vs.
indentation curve

0 .0.1 0.2 0.3 .0.4
Radial Coordinate r (mm)

Figure 3.13. Contact pressure distribution

0 0.2 0.4

Contact Radius (mm)

Figure 3.14. Contact force vs. contact radius










60 -

50 -

40 -

30 -

20 -

10 -


i, ,
- 4 mm plate //
- 3 mm plate /
- 2 mm plate

1 mm Plate /
i,' i
Eq. 3.13 // /
'// /

/ /

'/./ /

oo I W

I l I I
0 0.2 0.4
Contact Radius (mm)

Figure 3.15. Contact force vs. contact radius, comparison



0 0.01 0.02 0.03
Center Deflection (mm)

Figure 3.16. Contact force vs. center deflection

80 -
70 _____________
2 ---- Vertical Force (No Slip)
0o- Vertical Force (No Friction)
60 -



O 40
a /

S40 /

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016
Center Deflection (mm)
Center Deflection (mm)

Figure 3.17. Friction effect



4.1 Failure Criteria

The calculated stresses and strains from the axisymmetric

elements are based on the assumption that the plate is

transversely isotropic. However, the plate consists of

multiple plies where each ply is orthotropic and has different

fiber orientation. Hence, direct application of any existing

failure criteria which are based on stresses in the material

coordinates is not possible. An estimate of stresses in the

material coordinates has to be obtained from the axisymmetric


The existing failure criteria for composite materials can

be divided into two types. One is a interactive failure

criterion and the other is a set of independent failure

criteria for different failure modes. The interactive failure

criterion is basically an empirical technique, and

determination of the components in the polynomials of the

criterion requires difficult and expensive biaxial tests. The

most important criterion in this category is the Tsai-Wu

theory, which uses a general quadratic polynomial including


linear terms. This theory says that a failure surface in the

stress space exists in the form as

Fjoj + Fj = 1 ij = 1,.. 6 (4.1)

where Fi and F.. are second and fourth order strength tensors.

This criterion can not be widely used because of difficulties

in obtaining Fij. These criterion only predicts the onset of

failure and not the failure modes. However, if Fi and Fj can

be obtained correctly, this criterion may be satisfactory, and

it is simple to use [13,14].

Independent failure criteria, typically maximum stress

theory or maximum strain theory, can predict the onset of

failure and also failure modes. Distinguishing failure modes

is very important for the design of composite materials and

also for the progressive finite element program which must

modify its total stiffness matrix depending on the failure

modes. However, it lacks accounting for the interaction of

various components of stresses or strains in the failure


Good failure criteria must predict the principal failure

modes separately, such as fiber breakages, matrix cracking and

delamination, and each failure mode criterion should be

interactive such as the quadratic formula while they do not

require biaxial tests to determine the coefficients. The

quadratic failure criteria must take into account of the


compressive or tensile stresses. The work of Hashin, which

deals with the three-dimensional stresses fits this category

and is given in Ref. [13]. It leads to six independent

criteria, given in Equations (4.2) to (4.7). As input to the

failure criteria, the following material strength parameters

are needed.

X, ; fiber direction tensile strength

Xc ; fiber direction compressive strength

YT ; transverse tensile strength

Yc ; transverse compressive strength
S12 ; shear strength in the 1-2 plane

S1, ; shear strength in the 1-3 plane

S23 ; shear strength in the 2-3 plane

X, ; interlaminar (through-the-thickness) tensile stress

S1 ; interlaminar shear strength

Fiber failure criteria have two modes; compressive fiber

buckling mode and tensile fiber breakage mode.

For tensile fiber direction stress, ao > 0, Hashin's fiber

failure criterion is

S + 12 + 13 1 (4.2)
T S 12 S13

For the fiber buckling mode, oa < 0, the criterion is

S 1 (4.3)

There are also two types of failure modes for the matrix
cracking, depending upon o2 + 03, such that when (02 + 03) > 0,
the matrix cracking criterion is

(2 + (Ti2 2 ( 1(434)
23 + 23-0203 ) + 12 =1 (4.4)
YT S23S12 S13

and when (02 + 03) < 0, the matrix cracking criterion is

1 -1 (o02 +) + 2( 1-- (2+o 3)2
Yc[ 2S23 4S 23
T23 -02 03 12 3 )2
S23 S2 3 S13

There are two basic approaches to predict delamination
initiation [32]. One is the mechanics of materials approach,
which essentially compares the local state of stress in the
interply matrix layer where delamination occurs with relevant
strength parameters. The other approach is based on the
application of fracture mechanics principles and will not be


discussed here. However, once the delamination initiates, the

concept of a critical strain energy release rate, Gc, will be

applied to predict delamination growth.

In this study, delamination was predicted by the

mechanics of materials with a simple quadratic interactive

formula which depends upon 03, such that when 03 < 0, only

shear stresses induce delamination such as Mode II type, and

the criterion is

2 2
T13 +2 =2 1 (4.6)

when o3 > 0, it behaves like a combination of Mode I and

Mode II type crack initiation, and the criterion is

S)2 T13 + 23 =1 (4.7)
x, ) SI

Here XI and S, can be assumed as the tensile strength and

pure shear strength of the matrix. This assumption is

justified by the experimental observations that toughened

matrix materials apparently reduce the onset of delaminations.

Assume damage of a particular nature occurs at a location

given by a in the 00 layer. Then the same type and extent of

damage will occur in the 450 layer at (450 + a) and in the -450


layer at (-450 + a) and so on. That means in a quasi-isotropic

laminate the damage pattern will look like a stair-case in the

form of a helix. This phenomenon has been observed in impact

experiments also [15,33].

At a particular segment where the fiber direction is at

angle 0 to the x axis, the stress components in material

coordinates can be obtained as follows

01 cos2 sin28 2cos6sin8 0 0 0 o
02 sin28 cos26 -2cos0sin8 0 0 0
T12 -cosSsinO cosOsinO cos20-sin2 0 0 0 on
03 0 0 0 1 0 0 Ozz
T13 0 0 0 0 cos8 sin0 oa
3 0 0 0 0 -sinO cosO oez


because oae and Oez are zero for the axisymmetric analysis,

stresses can be easily obtained as in Equation (4.9).

o0 = cos2 o r + sin2 8 Oe

02 = sin2 08 a + cos2 8 ge

tz2 = cos 0 sin 8 (or Ge)
03 = zz

T13 = COS 0 orz

Tz3 = sin 8 orz


By combining Equations (4.2) through (4.7) and Equation
(4.9), the failure criteria for implementation into the finite
element program are obtained as follows.

Fiber failure

(i) cos28 or + sin28 ao > 0

cos2 o0 + sin2eo) O' cos) sinO (o-Or) 2 Cos0 rz2
XT S12 S13

( ii ) cos2 or + sin2e a < 0

cos28 r + sin28e

Matrix failure

(i) sin28 or + cos28 o + ozz > 0

1 (sin2 z sinr or zz-CS2 0 o)
cos9 sin9 (a0 -0a) 2 cos rz \2
S12 ( 13)

+ sin2 or + cos2 O + azz
YT 1

(ii) sin28 ar + cos2e o + azz < 0

-(sin2e rz-sin Ozz-cOs2e 09 ozz)

Scos9 sine (ae-o,) cose arz
S12 S 13

+ 4 (sin28 O +cos28 ae+6zz)2 = 1
4 53'


(i) ozz < 0


(ii) ozz > 0

orz)2 z)2
SI I = ) 1

Note that the delamination failure criteria are not

functions of 8. This fact can also be verified from the C-scan

results shown in Figure 1.5.

0 Ply

One Element

45 Ply

-45 Ply
-45 Ply

0 Ply
90 Ply

Figure 4.1 Partial failure patterns

Table 4.1. Ultimate strength of AS4/3501-6 lamina
unit (MPa)

XT ; 2172 Xc ; 1724

YT ; 53.8 Yc ; 205.5

S12 ; 120.7 S 3 ; 120.7

S23 ; 89.3 X ; 56.5

S, ; 89.3

* Courtesy of Hercules Aerospace Company


4.2 Implementing Delamination into the Finite Element Program

4.2.1 Introduction

Delamination initiates from where the shear stress is

high, such as above the mid-surface of the plate close to the

contact surface. Then delaminations start to expand radially,

and more delaminations appear in the bottom plies.

Delaminations can be simulated in the finite element program

by creating slip surfaces between elements as seen in Figure

4.2. New nodes are created, and rigid elements (spring

elements with very large spring constants K) are added between

two nodes across slip surfaces when the normal force is

compressive to prevent interpenetration of the nodes. Adding

new surfaces and consequently adding new nodes makes a

cumbersome procedure, which is a weakness of the finite

element method. The total stiffness matrix of the plate must

be recalculated after delaminations occur or grow. Also the

dimension of the total matrix is increased due to the

introduction of new nodes. If tensile stress exists at the tip

of delaminations, the spring element must be eliminated to

allow Mode I type (Opening Mode) delamination propagation.

Then, iteration of the program is necessary to assure that the

rigid elements are placed only where compressive stresses

exist. However, unless the contact force is high,

delaminations seemed to appear only at the upper portion of

the plate where through-the-thickness stresses oz are

compressive. The present finite element program monitors the


signs of all oa across the delamination surfaces while

implementing Hashin's failure criteria, and it was found that

all az's are compressive in the present study. It is safe to

assume that in this study the delamination growth is Mode II

type (Shear Mode) and to put rigid elements at all the nodes

across delamination surfaces. When delaminations become

extensive, especially beyond the mid-surface of the plate,

then the tip of the delaminations could have high tensile

stresses. This Mode I crack growth combined with Mode II type

will accelerate delamination growth faster than by Mode II

delamination alone.

4.2.2 Numerical Experiment

A simple numerical experiment was done to check the

effect of delaminations on the contact force vs. center

deflection curve and also to see the frictional effect across

the delamination surface. Assume that the experimental test

configuration shown in Figure 1.2 has pre-existing

delaminations of 12.7 mm diameter uniformly distributed

through the thickness. As 16 elements were employed through

the thickness, 15 delaminations were assumed in the finite

element program. An incremental point force was applied at the

center of the plate. Rigid elements were placed across all the

delamination surface nodes as shown in Figure 4.3.

Figure 4.4 shows contact force vs. center deflection

curve of the plate with pre-existing delaminations compared