Title Page
 Table of Contents
 List of Figures
 Analysis of interlaminar shear...
 Offset beam finite element for...
 On fracture behavior of a delaminated...
 Experimental study
 Conclusions and future research...
 Biographical sketch

Title: Dynamic delamination propagation in composite beams under impact
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00090194/00001
 Material Information
Title: Dynamic delamination propagation in composite beams under impact
Series Title: Dynamic delamination propagation in composite beams under impact
Physical Description: Book
Creator: Hu, Shoufeng,
 Record Information
Bibliographic ID: UF00090194
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 001577135
oclc - 22923152

Table of Contents
    Title Page
        Page i
        Page ii
    Table of Contents
        Page iii
        Page iv
    List of Figures
        Page v
        Page vi
        Page vii
        Page viii
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
    Analysis of interlaminar shear stress for an orthotropic beam due to static indentation
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
    Offset beam finite element for static and dynamic problems
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
    On fracture behavior of a delaminated composite beam subjected to impact
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
    Experimental study
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
    Conclusions and future research expectations
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
    Biographical sketch
        Page 139
        Page 140
        Page 141
Full Text








The author would like to express his sincere

appreciation to his advisor Dr. Bhavani V. Sankar. For his

most valuable assistance and constructive criticism, the

author is deeply indebted. Appreciation is extended to

Professors Beatty, Malvern, Sun, and Taylor for their

assistance and kindly serving on the author's doctoral

supervisory committee. The author would also like to thank

Mr. Young Kwon for his help in carrying out experiments.

In addition, the author would like to acknowledge the

sponsorship of the Ministry of Education of the People's

Republic of China and the Beijing Institute of Aeronautics

and Astronautics for his first year of study. Acknowledge-

ment is also due to the Department of Aerospace Engineering,

Mechanics, and Engineering Science at the University of

Florida for its financial aid during the subsequent years.

Without their support, the author's trip to the United

States and the continuous study would not have been


Finally, the author is grateful to his family, the

one he loves, and all his friends, whose encouragement, love

and friendship have always been the major motivation for the

completion of this work.



ACKNOWLEDGMENTS ...................................... ii

LIST OF FIGURES ...................................... v

ABSTRACT .......... ................................ vii


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

1.1 Background ................................... 2
1.2 Objective of This Study ....................... 4


2.1 Introduction ........... ....................... 6
2.2 Elasticity Solution ........................... 8
2.3 Numerical Results ............................ 16
2.4 Further Discussion of Stress Distributions .... 21
2.4.1 Contact behavior ......................... 18
2.4.2 Load-contact length relations ........... 22
2.4.3 Shear stress investigation ............... 23
2.4 Conclusions ................................... 29

DYNAMIC PROBLEMS .............................. 30

3.1 Introduction ................................ 30
3.2 Finite Element Formulation .................... 31
3.2.1 Method of weighted residuals --
a differential formulation ............... 36
3.2.2 Principle of total potential energy --
a variational formulation ................ 41
3.2.3 Lumped mass matrix ....................... 46
3.3 Performance of Offset Beam Finite Element ..... 48
3.3.1 Example 1: displacement and rotation of
tip of the beam ......................... 50
3.3.2 Example 2: strain energy release rate of
an end notched beam ...................... 51
3.3.3 Example 3: natural frequencies of
the beam ................................. 55


3.4 Newmark Method ................................ 58


4.1 Introdution ................................... 61
4.2 Discussion of Energy Release Rate ............ 63
4.3 J-Integral and Crack Closure Integral ........ 67
4.4 Rigid Element ................................. 72
4.5 Elastic Foundation Spring Model ............... 78
4.5.1 Principle of energy balance............... 79
4.5.2 Investigation of spring constants ....... 81
4.5.3 Dynamic crack propagation ................. 89
4.6 Buckling Effect .............................. 92
4.7 Conclusions and Discussion .................... 93

5 EXPERIMENTAL STUDY ............................ 97

5.1 Introduction .................................. 97
5.2 Apparatus and Specimen ..................... 97
5.3 Impact Test Procedure ........................ 100
5.4 Theoretical Prediction and Discussion ......... 103
5.4.1 Impact force history .................... 103
5.4.2 Dynamic crack propagation ................ 107
5.5 Concluding Remarks ............................ 109


6.1 Conclusions ............. .................. .. 113
6.2 Future Research Expectations .................. 115


METHOD .......... ............................. 118

SPRING MODEL ...... ........................... 120


HISTORY ........................................... 130

D.1 Effective Mass of a Simply Supported Beam
Under Impact ................................ 130
D.2 Determination of Impact Force History ......... 131

REFERENCES ........... ................................ 135

BIOGRAPHICAL SKETCH ................................. 139


Figure Page

2.1 A homogeneous and orthotropic beam
subjected to a uniformly distributed
load over a length 2a .......................... 9

2.2 Superposition of (a) half-plane solution,
(b) Timoshenko elastic beam solution, and
(c) quasi-elastic solution ................... 12

2.3 Superposition of the three cases for stresses
of a short beam due to three point bending
with distributed load and supports ............. 15

2.4 Normalized transverse shear stress r
distributions at x = 0.05, 6.25, and 12.0 mm ... 18

2.5 Normalized normal stress ax distributions
at x = 0. and 12.5 mm .......................... 19

2.6 Normalized normal stress a distributions
at x = 0. and 12.5 mm ......................... 20

2.7 Normalized maximum interlaminar shear
stresses for beams with different materials
under different loads .......................... 25

2.8 Locations of maximum shear stress from top
surface of the beam for different materials
under different loads ........................ 26

2.9 Comparison of normalized shear stress
distributions between uniform and elliptical
loadings ....................................... 28

3.1 Configuration of top and bottom beam
elements with nodes shifted .................... 32

3.2 Configuration of a sample cantilever beam ...... 50

3.3 Equivalent cases for calculating
displacement and strain energy release rate .... 53

4.1 Finite element nodes at crack tip .............. 70

4.2 Dimensions of the impact specimen used in
Ref. 6 ........................................... 71

4.3 Energy release rate by zero volume integral ..... 73

4.4 Energy release rate by crack closure integral ... 74

4.5 Transverse displacement at left crack tip ....... 75

4.6 Comparison of energy release rates by J-
integral, continuous spring model, and discrete
spring model for Mode I fracture ................ 83

4.7 Energy release rate by spring model ............. 87

4.8 Energy release rates by three different
approaches .......... ............................ 88

4.9 Experimental and finite element results of
left crack tip position: crack extension
versus time ..................................... 91

4.10 Spring model in finite element analysis for
a delaminated plate ............................. 96

5.1 Pendulum impact test facility scheme ............ 98

5.2 Dimensions of the specimen for pendulum impact .. 101

5.3 Dynamic responses for specimens with and
without crack extension during impact: impact
force versus time ............................. 105

5.4 Left and right crack extensions versus time ..... 110

B.1 Spring model for Mode I fracture (symmetric) .... 122

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



Shoufeng Hu

May, 1990

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

The major objective of this study is to analyze

interlaminar crack growth of composite beams with embedded

delaminations under impact loading.

First, a superposition technique is proposed for

obtaining an analytical solution for the shear stress field

in a short composite beam under transverse loading, for

which the shear stress plays a significant role in the

interlaminar failure. The results calculated are accurate

and have shown a high shear stress concentration. The

delamination is most likely to initiate at the location

where this high shear stress occurs.

A node-shifting technique is used to solve two-

dimensional fracture problems using a one-dimensional beam

element. The technique has been successful in calculating

the energy release rate of a beam with an embedded


delamination and achieved savings in computer time and

storage. This is meaningful, especially for the dynamic

problems with a large number of degrees of freedom.

An elastic foundation spring model is used to compute

the dynamic energy release rate. In the model a delaminated

beam is replaced by two separate beams joined together with

a distributed spring at the uncracked part of the original

beam. The energy release rate of the crack is equivalent to

the potential energy per unit area of the spring at the

crack tip. This principle also can be extended to the finite

element analysis. In this study the spring model has been

extensively investigated and compared with other models.

With an appropriate choice of the spring constants,

identical results for the energy release rates have been

obtained by J-integral, crack closure integral and the

spring model.

Finally pendulum impact tests have been performed to

investigate the validity of the theoretical model and the

numerical analysis, as well as to compare pendulum impact

with projectile impact. A satisfactory agreement between

experimental and theoretical results has been obtained. Some

distinctive phenomena of pendulum and projectile impacts are

also discussed.




Improved specific strength and stiffness properties

have made fiber-reinforced composites very attractive for

use as advanced structural materials in the aerospace

industry. Composite materials have been receiving increasing

attention not only because they are ideal for the weight

sensitive structures such as aircraft, rockets and

spacecraft, but because their anisotropic property can be

utilized to achieve enhanced performances. The forward-

swept wing is an example which expanded horizons in the

application of new materials.

However, unlike most metals, the damage tolerance of

composite materials is limited and the performance in

failure is different, since the composite materials commonly

used for aircraft structures (e.g. graphite/epoxy) are

brittle and have no plastic stress relief capability. For

example, even in normal use an aircraft may be subjected to

a low energy impact from dropped service equipment, which

causes no visible surface damage. For conventional

structures such minor incidents may produce some slight

surface dents which will not seriously reduce the strength

of the structure. But for composites, these impacts may

result in interlaminar cracks and these cracks may propagate

under certain conditions such as fatigue loads, impact loads

or change of temperature. Hence they may cause buckling,

partial structural failure and finally failure of the entire

structure. Since interlaminar crack is one of the most

commonly observed failure modes in composite structures, the

study of it is of great necessity and importance in fracture

analysis and failure prediction for composite structures.

1.1 Background

Efforts have been made in the past few years to

investigate the crack problems of composite materials.

Interlaminar and intralaminar fractures by tension,

compression and transverse loading have already caught

researchers' attention. Mode I type fracture, believed to be

a matrix dominated fracture, occurs in cases like central-

notched plates in tension. Mode II type fracture, controlled

by the fiber-matrix interface or the interface between two

layers, occurs when plates or beams, e.g., subjected to a

bending, twisting or transverse impact load.

Since the interlaminar crack propagation is dominated

by Mode II type fracture in composite plates under an impact

load, studies have focused on Mode II type static and

dynamic fracture analysis. S. S. Wang [1], A. S. D. Wang and

Crossman [2] initiated research work about ten years ago on

the mechanism of delamination, the initiation and growth of

transverse crack and edge delamination in composite

laminates under static loads. Taking the end notched flexure

(ENF) specimen as a candidate for measuring the strain

energy release rate, Carlsson et al. [3] examined the

influence of the interlaminar shear deformation, the

friction between the crack surfaces, and the geometric

configuration on the strain energy release rate. Lagace [4]

and Takeda et al. [5] have worked on the interlaminar crack

nucleation of composite plates under impact analytically and

experimentally. Their independent efforts were basically

directed toward studying the response of composites to

impact, analysis of the impact event, and the damage induced

by impact.

Concentrating on experimental work, Grady and Sun [6]

first investigated the phenomenon of the propagation of an

existing crack under an impact load. In their work, the low

energy impact tests were performed by impacting the

specimens with a low mass and high speed projectile, and

recording the impact response and subsequent crack

propagation by high speed photography. From the photographic

data, the contact duration and the dynamic response of the

position of the specimen were measured. The velocity of the

crack was evaluated. Finally, the time and position

dependent nature of the crack growth velocity and its

variation with impact conditions were investigated.

1.2 Objective and Outline of This Study

Since it is not always possible to examine the

cracks, voids, or any kind of damage around the entire

aircraft structure after each flight, it is necessary to

allow for the existence of cracks and to ensure that these

cracks will not propagate and consequently cause any

structural failure during the service period. Therefore, in

design, the existence of cracks should be taken into

consideration. The damage tolerance of a structure with

existing cracks should be an integral part of the design. It

is the primary purpose of this study to investigate and

understand the mechanism of the existing interlaminar crack

propagation under an impact load, and finally to prevent it.

This study consists of six chapters. In chapter 2, a

superposition technique is proposed to analyze the stress

fields in a short laminated composite beam under a static

load, for which high stress is the major cause of

interlaminar failure. Chapters 3 and 4, in turn, are devoted

to finite element formulation and investigation of fracture

behavior for a delaminated beam. An effective and accurate

offset beam finite element is developed. It is used, in

conjunction with the spring model covered in chapter 4, to

solve for the energy release rate and crack growth of a

delaminated composite beam under a transverse impact load.

Chapter 5 is focused on experimental study. Pendulum impact

tests are performed to verify some aspects of the


theoretical models. Chapter 6 includes the summary of the

present study, concluding remarks, and suggestions for

future studies.



2.1 Introduction

Composite laminates in transverse impact are highly

susceptible to delamination due to high shear stress near

the contact area. Delamination is one of most commonly

observed failure modes in composite laminates under impact

loading, especially if the laminates are small in dimension.

It occurs whenever the highest interlaminar shear stress

reaches or exceeds the strength of the laminar interface,

which is considerably weaker than the fiber/matrix

interface. However, it has been reported that the stress

distributions in orthotropic beams under a transverse load

are significantly different from beam theory results [7].

The contact behavior between a rigid indentor (impactor) and

an orthotropic beam may also make the displacement and

stress responses very complex. Therefore estimation of

actual contact area, of contact stresses therein, and the

detailed stress field in orthotropic beams is essential for

a comprehensive understanding of failure mechanisms, and a

meaningful interpretation of experimental results.

Further, the location of the maximum shear stress

will most likely be the position where the interlaminar

crack initiates. The stresses in the vicinity of impact can

be approximated by corresponding values of the static

contact problems, since the elastic wave effect is presumed

to be small due to the dimension. Hence the study of a beam

under static loading can be led to obtain the information

regarding the nature of interlaminar shear stress


A simple, yet accurate procedure is developed in the

present study for obtaining an analytical solution for

stresses in an orthotropic beam under three-point bending,

which is not directly solvable by either classical beam

theory or the available elasticity methods. Attention is

focused on contact behavior, contact stress distribution and

shear stress field in the vicinity of contact for different

materials. A superposition method developed in the present

study combines an orthotropic half-plane solution, a

Timoshenko elasticity solution for beams and a quasi-

elasticity solution. A good agreement among the results of

Whitney [7], finite element analysis and this study has been


With the aid of the contact theory, the stress

distributions, especially the interlaminar shear stress, of

an indented beam are calculated by the superposition method.

The maximum shear stresses for different load distributions

and contact areas are examined. The effect of material

properties on the shear stress is discussed.

2.2 Elasticity Solution

In this section a homogeneous orthotropic beam

subjected to a load P uniformly distributed over a length 2a

(Fig. 2.la) is introduced. The assumptions are that the

material principal axes are parallel to the coordinate axes,

the beam length 2L is much greater than the thickness h, and

two supports at the ends are distributed shear stresses

along the beam thickness. The stresses in this simply

supported beam are computed by superposing the stresses in

an orthotropic half-plane, an elastic beam suggested by

Timoshenko and a quasi-elasticity beam. In the following

paragraphs the notations a,, a2 and a3 are used to denote

aXX' yy and rxy respectively. The beam width is assumed to

be unity.1

The stresses in an orthotropic half-plane due to a

unit normal concentrated force at the origin are given by

Leknitskii [8] as

arr = -(pa + p2)Jf 1122CosO/irrL(0)

a08 = arO = 0 (2.1)

iIn order to compare the results with that from Whitney
[7] we use the plane stress assumption. In reality, for the
investigation of interlaminar shear stress only plane strain
assumption is meaningful. Whitney [7] indicated that for
these commonly used material properties the differences
between the compliance coefficients of plane stress and
plane stain are negligible. Therefore the results of the
plane stress assumption will not significantly differ from
that of the plane strain assumption.

i ^2a




1: h

-- Cross Sections
for Normal Stresses

.......... Cross Sections
for Shear stress

Unit: mm

Figure 2.1 A homogeneous and orthotropic beam subject to a uniformly
distributed load over a length 2a






where L(0) = fi1sin4 + (2812+Pg6 )sin2 cos2 + p22cos4 ,

and pi, A2 are the roots of the characteristic equation

p1114 (2812+P66)p2 + 822 = 0 (2.2)

For plane stress pij = Sij and for plane strain ij =

Sij -S13Sj3/S33. The compliance coefficients Sij are related

to the material elastic constants as S I = 1/El, S22 = 1/E2,

S33 = 1/E3 S12 = -v12/E S13 = -v13/E1, S23 = -v23/E2 S66

= 1/G12, and Sij = Sj The cartesian components of

stresses in Eq. 2.1 are

ax = gi (x,y) = ar sin20

ayy = 92 (x,y) = arCOS28

Txy = 93 (x,y) = a,,sinOcoso

The stresses due to the uniformly distributed normal

force P/2a on the edge of the half-plane can be calculated

from the point force solution as

a (x,y) = (P/2a) gi (x-C,y)dC, i = 1,2,3 (2.3)

where the superscript h stands for half-plane.

Since the bottom surface of the original beam at y =

h is traction free, we tend to apply a normal and a shear

traction equal and opposite to ac (x,h) and ah (x,h)

respectively to the bottom of the beam, and superpose the

resulting stresses to the half-plane solution (Fig. 2.2).

The length of the beam is assumed to be infinity. The

reactions are located at x = L. In the following paragraphs

a Timoshenko elasticity solution for beams for stresses due

to the residual traction -a (x,h) and a quasi-elasticity

solution due to -a (x,h) are derived.

Consider a simply supported orthotropic beam

subjected to a distributed transverse load q(x) = -ah (x,h)

on the bottom surface y = h (Fig. 2.2). The stresses in the

beam can be obtained as

ai (x,y) = ai (classical beam solution) + {corrections) (2.4)

where the superscript bl represents the Timoshenko

elasticity solution for beams. Then an approximation for the

correction terms can be obtained as explained below.

The exact solution for the stresses in a simply

supported isotropic beam subjected to a uniformly

distributed load, say q0, along the bottom surface y = h is

given by Timoshenko and Goodier [9]. Using the same

procedure, the stresses in an orthotropic beam can be

determined as


I -1

(b) (c)

Figure 2.2 Superposition of (a) half-plane solution, (b) Timoshenko elastic
beam solution, (c) quasi-elastic solution


oax = 12M(x)(y 0.5)/h2 + {aq0 (4y3 6y2 + 2.4y 0.2)}

ayy = {q (-2y3 + 3y2))

Xy = 6Q(x)y(l y)/h (2.5)

where M(x) is the bending moment, Q(x) is the shear force,

y = y/h, a = (2#12+P66)/2pi,1 and terms inside {} represent

the corrections to the classical beam theory solution. If,

as in the present case, the load is not uniform but does not

vary significantly, the correction terms can be approximated

by substituting q(x) for q0 in the terms inside () in Eq.

2.5. Such an approximation for isotropic beams is proposed

by Timoshenko and Goodier [9]. The bending moment and the

shear force can be found as

M(x) = P(L x)/2 Jq() (Q x)dC (2.6)
Q(x) = -P/2 + Lq(()d (2.7)

Now consider the case of shear traction t(x) =

-a (x,h) on the bottom surface y = h. Because of symmetry

the shear forces are self-equilibrating and there will be no

support reactions. A solution for stresses that satisfies

the equilibrium equations and boundary conditions on the top

and bottom is derived as

axx (x,y) = ao2 = 2t(x) (1 3y)/h (2.8)

ay (x,y) = ab2 = t'(x)hy2 ( y) (2.9)

xy (x,y) = ab2 = t(x)y(3y 2) (2.10)

where tl (x) = t(Q)d(, t'(x) = dt/dx and the superscript

b2 represents the quasi-elasticity solution.

This solution does not satisfy the compatibility

equations exactly but does so in an approximate sense if

t(x) is a slowly varying function of x, i.e., t'(x) = 0. It

may be noticed that Eq. 2.8 is the exact elasticity solution

for an infinitely long beam, whereas Eqns. 2.9 and 2.10 are

obtained by integrating the differential equations of


In conclusion, we have derived expressions for ah,

a1 and a2 which superposed together represent the stress

fields in an orthotropic beam subjected to a transverse

load P distributed uniformly over a length 2a, i.e.,

ab = a bl + a\2 (2.11)

The detailed expressions for ah a o and aC2 are

given in appendix A.

Referring to Fig. 2.3 the stresses under three-point

bending then can be obtained by superposing stress again as



Figure 2.3 Superposition of the three cases for stresses of a short beam

due to three point bending with distributed load and supports


| ...

ai (x,y) = Lo (x,y) + 0.5aL (-x-d, h-y)

+ 0.5a4 (d-x, h-y) (2.12)

where 2d is the span between two distributed supports and it

is also the effective beam length. The superscript L is for

the stresses obtained from Eq. 2.11 with an assumption of a

long beam. By this superposition, the stress distributions

are obtained for short beams, which are more susceptible to

interlaminar shear failure.

Finally all stresses are non-dimensionalized as

follows, in order to be compared with the stresses obtained

by whitney [7]:

Bh2 2Ba Bh
axx 2Pd-xx yy p ayy xy prxy

2.3 Numerical Results

Whitney [7] performed a two-dimensional Fourier

series elasticity analysis for an orthotropic beam under

three-point (and four-point) bending, where the concentrated

load and supports are represented by a uniform normal

traction distributed over a small length of the beam. In his

paper, a set of numerical results has been published, which

can be used to compare with and examine the accuracy of the

results of this study.

ASTM Standard D-2344 [10] provides the standard for

the interlaminar shear test of a short composite beam. In

Ref. 7 some of the results are obtained with the

configuration of this standard shown in Fig. 2.1b. The

material constants used (graphite/epoxy) are as follows:

EI = 145 GPa E2 = E3 = 10.3 GPa

V12 = 13 = 0.3 v23 = 0.55 G12 = 5.52 GPa

A finite element analysis (Eight-Noded Isoparametric

Two-dimensional Quadrilateral Parabolic Element) is

performed to calculate the stress distributions in order to

verify the results of Whitney's study and the present study.

Figures 2.4 and 2.5 show a good general agreement

among Whitney's Fourier series solution, the closed-form

solution of the present study and the finite element

numerical solution. The ayy predicted by the present study

differ from Whitney's results. However they agree very well

with that of the finite element solution (Fig. 2.6).

2.4 Further Discussion of Stress Distributions

The investigation of the contact between a rigid

cylinder as an indentor and an orthotropic beam plays an

important role in evaluating contact force, contact area and

contact stress distribution over an orthotropic beam. This

problem has been previously solved by Sun and Sankar [11].

2.4 -

2.2 -- This Study

2 A Ref.7


1.4 x=0.05 mm

S 1.2 x=12.0 mm

1 \ x=6.25 mm

0.8 -

o0. -

0.4 -

0.2 -

0 2 4

Y-axis (mm)

Figure 2.4 Normalized transverse shear stress Tau(xy) distributions at

x = 0.05, 6.25, and 12.0 mm



-1 -

2 -2 -
0X0 11
E -3 xmm x=12.5mm

D -4

c -5 +
-6 -5

cc -8

z -9 This Study

-10 -

-11 -

0 2 4 6
Y- axis (mm)

Figure 2.5 Normalized normal stress Sigma(x) distributions at

x = 0. and 12.5 mm

0. -
0.9 Ref.7

-- This Study
D 0.8

~0 o.

Z x=O. mm x=12.5 mm

Z 0.2

0- II

0 2 4 8
Y- axis (mm)

Figure 2.6 Normalized normal stress Sigma(y) distributions at

x = 0. and 12.5 mm

In this section, we first briefly describe the contact

behavior and then outline its applications to the present


2.4.1 Contact behavior

For an indentor with a small radius in comparison

with the thickness of the beam, i.e. a/h is small (here a is

half the contact length), the normal contact stress

distribution can be obtained by Hertzian orthotropic half-

plane solution [12] as

p(x) = (2P/?a)Jl (x/a)2 (2.13)

If the radius of the indentor is very small compared

to the average radius of curvature of the deflected beam,

the contact behavior will be essentially the same as that of

the orthotropic half-plane.

For a relatively large radius indentor or a thin

beam, a numerical evaluation process has been provided by

Sankar [13]. It has been observed that a large ratio of a/h

will result in a stress distribution very different from the

Hertzian solution.

In his recent work, Sankar [14] developed a closed

form formula to modify Hertzian orthotropic half-plane

solution as a function of a/h,

p(x) = (2P/7a)/1 x2 [1 f(i 4x2)] (2.14)

where x = x/a, p = 8.75B(a/2d)4 and

B = r/16(E2/E )(l/p1+p2)(2d/h)3 .

It can be directly observed how these different

constants influence the contact stress distribution. If f is

very small, the contact stress distribution is close to an

elliptic shape, i.e. Eq. 2.14 degenerates back to Eq. 2.13.

If the factor p is relatively large, the contact stress

distribution will take the shape of a butterfly [12].

2.4.2 Load-contact length relations

Similar to the situation discussed above, for a

relatively small radius indentor (or a thick beam), the

load-contact length relation can be derived to be the same

as that for an orthotropic half-plane [12],

P = KCa2 (2.15)

where KC = 1/2KhR, Kh = (p1+p2)/7rE2, i1 and ,2 are the roots

of Eq. 2.2 and R is the radius of indentor [8, 12]. For a

relatively large indentor or a thin beam the curvature of

the flexed beam will be comparable to that of the indentor

and hence will have a significant effect on the contact

area. Therefore, it can be considered as the contact between

two curved bodies, and then Eq. 2.15 will be modified as


P = K, (1 R/Rb)a2 (2.16)

where Rb is the average beam radius of curvature at the

center given as Rb = 4EI/Pd.

For a typical contact problem with the contact force

given, the contact length and then the contact stress

distribution can be obtained by Eqns. 2.13 and 2.15 for

small ratio of a/h, or by Eq. 2.14 and 2.16 for a large

ratio of a/h. The contact stress can be replaced by the sum

of a number of uniformly distributed loads over varying

lengths symmetrical with respect to the center. The stress

field due to each uniformly distributed load can be obtained

by the superposition method explained in section 2.2. Then

the stress field due to the contact stress expressed as Eq.

2.13 or 2.14 can be obtained with the sum of each stress

field due to the uniformly distributed load.

2.4.3 Shear stress investigation

Three kinds of materials used in the this study [15]

are in Table 2.1.

Table 2.1 Engineering constants for typical unidirectional

Type Material El E2 V12 G12

T300 Graphite
/5208 /Epoxy 181 10.3 0.28 7.17

Scotchply Glass
/1002 /Epoxy 38.6 8.27 0.26 4.14

B(4) Boron
/5505 /Epoxy 204 18.5 0.23 5.59

ASTM standard 2344 is used for the specimen

configuration. All the stresses are normalized. It has been

found that the maximum shear stress always occurs at x = a,

the end of the contact region irrespective of the contact

stress distribution. However, contact force, contact area,

radius of indentor, thickness and span of the beam, and

material properties all are the influential factors for the

contact stress and thereby the magnitude and the location of

the maximum shear stress. For the standard beam

configuration (Fig. 2.1b) and radius of indentor (3.175 mm),

the different longitudinal and lateral Young's moduli El and

Ez result in a different contact length and contact stress

distribution, and thereby a different magnitude and location

of the maximum shear stress, see Figs. 2.7 and 2.8 (the

contact forces are 50 N, 100 N, 200 N, 500 N, and 1000 N

respectively). The magnitude and location of the maximum

shear stress are very complex functions of material

constants E, and E2. It can not be directly observed if a



6 + Glass/Epoxy

S5 Boron/Epoxy

E 4

3 -

+ +
0 +

,.- i ,i-- i ,-- i i -- ,- i -- ,- i -- ,- ,- -- -- -- i i-- --
0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
Contact Length a (mm)

Figure 2.7 Normalized maximum interlaminar shear stresses for beams with

different materials under different loads



1.3 0 Graphite/Epoxy

1.2 -
+ Glass/Epoxy
S 11 -
1 Boron/Epoxy

S0.9 -
C 0.8 -
E 0.7 -

M 0.6 -

S 0.5 -


0.3 +

0.2 +
+ M
0 .1 i I i I I I I I
0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
Contact Length a (mm)

Figure 2.8 Locations of maximum shear stress from top surface of

the beam for different materials under different loads

larger E2 or a smaller El will produce a smaller maximum

shear stress.

For a given contact stress distribution, the trends

clearly indicate that a larger contact length produces a

smaller maximum shear stress, which is located further away

from the top surface of the beam. Furthermore different

shapes of the contact stress distribution will alter the

magnitude and the location of the maximum shear stress. For

the same contact force (50 N) and contact length (0.4206

mm), Fig. 2.9 illustrates that an elliptically distributed

contact stress produces a smaller maximum shear stress than

that produced by a uniformly distributed contact stress. The

comparison is for the shear stress distributions between the

top surface (y=0 mm) and the position y=1.0 mm. Below that

position the difference is negligible. In contrast, a

butterfly shape contact stress may produce a larger maximum

shear stress. However, the butterfly shape contact stress

distribution only happens for a very large indentor on a

very thin and long beam. For the dimensions shown in Fig.

2.1b, the radius of indentor has to be as large as the

length of the beam to result in a butterfly shape contact

stress. Under this condition the shear stress loses its

significance in failure as well as the significance for



0 Uniform Loading

A Elliptical Loading

e Beam Material: Graphite/Epoxy
Load: 50 N
S 5 -Contact Length: 0.4206 mm

(D 4


z 2

0 0.2 0.4 0.6 0.8
Y axis (mm)

Figure 2.9 Comparison of normalized shear stress distributions

between uniform and elliptical loadings

2.4 Conclusions

The superposition method developed in this chapter

can be used for beams with various dimensions, loading

conditions and material properties. The whole procedure

provides a useful tool to evaluate the stress distributions,

as well as an access to understanding how these factors

alter the magnitude and location of the maximum shear


An important feature observed is that the shear

stress distribution is significantly skewed in the vicinity

of the contact area. In some cases the maximum shear

stresses could be more than ten times higher than that

obtained by classical beam theory.

With regard to the influential factors of the contact

between a rigid cylinder and an orthotropic beam, the

analysis presented here shows, in Figs. 2.7 and 2.8, that

based upon the same load, beam configuration and radius of

indentor, the position and the value of the maximum shear

stress can be significantly different for different

materials. In addition, not only load, contact area and

ratio a/h, but the shape of the contact stress distribution,

will influence the magnitude and location of the maximum

shear stress. Hence the assumption of uniform distribution

of the contact stress, which is commonly used, usually

overestimates the maximum shear stress.



3.1 Introduction

A one-dimensional beam element with nodes offset to

one side, either top or bottom, is introduced. Timoshenko

beam theory [16] is used to derive the relevant equations.

This shear deformable and node-shifted beam element provides

a simple and effective method for static and dynamic

fracture analyses of the laminated beams with interlaminar

cracks. The beam is assumed to be homogeneous, orthotropic

and linearly elastic. The material principal directions of

the beam coincide with the x-y-z axes. Shifting the nodes

from middle plane of the beam to the top or bottom surface

will introduce an extra shear force and a corresponding

axial displacement. The elements with nodes at the top and

those with nodes at the bottom join together to form a

complete beam. The delamination, which runs across the

width, is located between top and bottom elements where the

nodes are separated.

In comparison with the two-dimensional element model,

which is commonly used for fracture problems, the offset

beam element can greatly reduce the computing time and

computer storage, which is especially significant for

dynamic problems. This beam element may not be able to

calculate the stress intensity factor K, which characterizes

the crack tip stress field, with satisfactory accuracy. But

for computing the strain energy release rate or the J-

integral, which is related to the strain energy as a whole,

this beam element is powerful, and the results have shown a

good agreement with two-dimensional finite element results.

3.2 Finite Element Formulation

Based on the assumption of Timoshenko Beam Theory

which includes the shear deformation, we define i(x) as the

rotation about y-axis, and w(x) as the displacement of the

beam in z-direction (Fig. 3.1). Both of them are functions

of x only. Therefore,

ux = z#(x) (3.1)

u, = w(x) (3.2)

These two equations imply that planes which are

normal to the beam axis in the undeformed state remain plane

in the deformed state, and the vertical displacement of all

points on any normal cross-sectional plane is the same,

i.e., there is no thickness stretch. These are two

fundamental assumptions for beam theory.


X wu w wIu u

T- T

Figure 3.1 Configuration of top and bottom beam elements with
nodes shifted

Using the strain-displacement and stress-strain

relations, we readily obtain

S= Ez dx (3.3)


rxz = IG(O + d-) (3.4)

where E is the longitudinal Young's modulus and G is the

shear modulus in the x-z plane, which are equivalent to El

and G13 in composites.

In Eq. 3.3 the term v(ao+ao) is neglected and in Eq.

3.4 a shear coefficient K is introduced to adjust the

present beam theory to agree with the predictions of the

three-dimensional theory [16]. Cowper [17] provided a

detailed discussion of the determination of K. For a beam

with a rectangular cross-section and a Poisson's ratio of

0.3, K = 0.850 is suggested [16].

Upon substitution of Eqns. 3.3 and 3.4 into the

expressions for bending moment and transverse shear force we


M = [ZadA = EI- (3.5)
J dA = dxGA() (3.6)

Q = rzdA = xGA(1O + dw (3.6)
A dx

where I is the moment of inertia and A is the cross-section


Now consider the equilibrium equations for a

Timoshenko beam with longitudinal shear traction t(x) at the

top or bottom surface of the beam, which may respectively be

called bottom and top beam. Meanwhile, Eq. 3.1 becomes

ux = u (x) + z (x)

= u(x) + (z h/2)0 (3.1a)

where u0 (x) and u(x) are the displacements in longitudinal

direction respectively at the center and at the top (or

bottom) surface of the beam where shear traction is applied,

and h is the thickness of the whole beam and each of the top

or bottom beam shears the half of it. The positive sign is

taken for the bottom beam with the shear traction at its

top. Equation 3.5 is found to remain the same.

The first equation, which describes the equilibrium

in z-direction, is

dx + p(x) = 0

and it can be written in terms of displacements as

d-i GA(- + )] + p = 0 (3.7)

where p(x) is the external distributed load in transverse


The second equation, which corresponds to rotational

equilibrium, has different signs for top and bottom beam,

dM h
dx -Q (-t) = 0

The third equation is the equilibrium in x-direction

dx + t = 0

where N is the normal force in longitudinal axis of the beam


duo d h
N = EA =- EA(u + )
dx dx 2

Then these two equilibrium equations in terms of

displacements become

d (EIdw h d dG h
dx(EIdx iGA(dx + ) 2 dx EAdx(u + 2 = 0



dx EAd(u )] + t = 0 (3.9)

3.2.1 Method of weighted residuals-a differential

The method of weighted residuals is a technique for

obtaining approximate solutions to linear and nonlinear

partial differential equations. Applying the method of

weighted residuals involves basically two steps. The first

step is to assume the general functional behavior of the

dependent field variable in some way so as to approximately

satisfy the given differential equations, which are Eqns.

3.7, 3.8 and 3.9 in the present case, and boundary

conditions. Substitution of this approximation into the

original differential equations and boundary conditions then

results in some error called a residual. This residual is

required to vanish in some average sense over the entire

solution domain.1

The second step is to solve the equations resulting

from the first step and thereby specialize the general

functional form to a particular function, which then becomes

the approximate solution sought.

1In other words, the integral of the weighted residuals
over the entire solution domain is zero. The weighting
function can also be understood as a virtual displacement.
Then the method of weighted residuals is equivalent to the
virtual work principle.

To be more specific, we now consider the problem of
concern. By applying the method of weighted residuals, we

f d dw
vdx{nGA(dx + 0)] + p)dx = 0

d+d d -x dw hd du
0{ d(4EI ) iGA( + ) dx(EA d)}dx = 0

L d d h
f{dx[EAdx((U + ) + t}dx = 0

where v,o,f are weighted functions, or called virtual
displacements, corresponding to the vertical, rotational and
longitudinal directions respectively, and the terms inside
the parentheses are the residuals.
Integrating by parts, we obtain

GA( d + 7)dx = :pvdx + vxGA( + #)j

{4Edi -- + GAO( + ) + E du dx
od dx dax 2 dx
L h t
o 2 o

LEA(du df + h df)dx = ftdx + fEA(du + h d)L
o adx dx 2 dx dx fE dx 2 dx 0

The unknown displacements can be interpolated as

0 = DoiMi

and u = ui Ki ,


where the interpolation functions Ni, Mi and Ki are

interpreted as shape functions, and the unknown parameters

wi, i, and ui are defined as nodal displacements. In
particular, let these weighted functions be equal to one

term of these assumed shape functions, i.e., v = Nj, 0 = Mj

and f = Kj. Then the equilibrium equations for one-

dimensional offset beam finite element in matrix form are

I dNi dNj r GdNj
JGAAd~x dx j GAN Midx

d S dx d
iGAd i dx (4EdMidg + GAMiMj )dx EA dx dM-dx
fLdx f dx dx 2 dxd

fpNi dx


I dKi dK
L dx dx

J+ h dMidKdx
EL dx dx

Q(xL)Nj (xL) Q(xo)Nj (xo)

[M(xL)MJ (xL) M(xo)Mj (xo)
h h
[F(x )Mj (xL)2 -F(xo)MJ (xo)2]

F(XL)KJ (XL) F(Xo)Kj (Xo)

where L is the element length, xL and x0 are the nodal

coordinates of the right and left ends of the element, and

w = wi N ,

Q, M, and N are the forces corresponding to the

displacements w, 0, and u. Then the stiffness matrix is

obtained as

J.dNi dNj r dN
GAGGdNAdx GAa Midx
dx d d dx

dNi dMi dM h dKidMj
LGAaX M dx (4EI +GAMi )dx EA d dx dx

EA d dx EAdKdK dx
L 2 dx dx x dx

From Eqns. 3.7-3.9 one realizes that the function of

w is one order higher in differential form than the

functions of 4 and u. Therefore if one chooses shape

functions Mi and Ki as linear functions, then Ni, the shape

function for w, needs to be quadratic to avoid the shear

locking effect [18]. Hence, for each element there will be

three nodes. The nodes at two ends contain all the

displacement variables w, i, and u but the middle node

contains w only.

The displacements are interpolated as

w(x) = wiNi + w2N2 + w3N3

P(x) = 01MI + 03M3
u(x) = uK, + u3K3

Then the shape functions and their derivatives are

obtain as follows:

x x2 X
N = (1 2 (1 x) = 2 32 +1
L L L2 L

N2 = 4X(i ) =-4+ 42
L L' L L

N3 = (2 -1) = 22 -

dN1 4x 3
dx L2 L

dN2 x 4
dx -L2 L

dN3 x 1
dx L2 L

M, = KI = 1 -

M3 = K3 =

dMi dK= 1
dx dx L

dM3 = dK3_ 1
dx dx L (312)

Substitution of above shape functions and their

derivatives into Eq. 3.11 results in the general expression

of element stiffness matrix as

7 GGA 5 8 xGA 1 GA 1G
3 L 6 3 L 3 L 6GA

4EI L G EAh 2&GA 1 GA -4EI L -EAh
+,cGA + ,cGA ,cGA +GA +
L 3 2L 3 6 L 6 2L

L 2L L

16 xGA 8 KGA 2
3 L 3 L 3GA

7 KGA 5G
3 L 6

symmetric 4EI L EAh
+-mGA +
L 3 2L



the corresponding displacements, in turn, are w 1 ul,

W2 W3, 3 3 and u3

3.2.2 Principle of total potential energy-a variational

A continuum problem usually has different but

equivalent formulations-a differential formulation and a

variational formulation. In the differential formulation, as

we derived above, the problem is to integrate a system of

differential equations subject to given boundary conditions.

In the classical variational formulation the problem is to

find the unknown function (or functions) that extremizes

(usually minimizes) a functional (or functionals) subject to

the same boundary conditions. The two formulations are

equivalent because the functions that satisfy the

differential equations and their boundary conditions also

extremize (or make stationary) the functionals.

The classical variational formulation of a continuum

problem often has advantages over the differential

formulation from the standpoint of obtaining an approximate

solution. For example, the functional often actually

represents some physical quantity in the problem, like total

potential energy. Also, the functional contains derivatives

of order lower than that of the differential operator, and

consequently an approximate solution can be sought in a

large class of functions.

Consider the following energy conservation equation,

n(0) = W(4) U(O) T(O) (3.14)

where H, W, U, and T are total potential energy, potential

energy of external forces, strain energy, and kinetic energy

respectively, and 4 is the generalized displacement.

To extremize H with respect to 4 we require that

6 = a 6, = 0 (3.15)
i=1 aVi

where n is the total number of discrete values of 4 assigned

to the solution domain.

Since 64 's are independent (physically they are the

virtual displacements) and will not be zero, Eq. 3.15 then


0 i = 1,2,3,*--,n (3.16)

Upon the substitution of Eq. 3.14 into Eq. 3.16 we


8W_ d F aT + a 317)
+ (3.17)
a4i dt a 41 aji

This formula can also be recognized as Lagrange

equation. Moreover, the strain energy and the kinetic energy

can be written in the form

U = --Eki j ij (3.18)

T = 2 i (3.19)

where kij and mij are the components of the stiffness and

mass matrices.

Substitution of Eqns. 3.18 and 3.19 into 3.17 results

in the expressions for the stiffness and mass matrices as

dt ) = Xmil"

au k j
av j



Expressions for U and T in terms of nodal

displacements can be determined as

:Lh / 2 1
U = f (EIe2 + nGA-y )dZdx
S J-h / 2

T i pB(u2 + u2)dzdx
0 J -h/2



With the aid of strain-displacement relations

x ux
eX ax

Ux auz
'Yxz az + ax

as well as Eqns 3.la and 3.2, we obtain

U = 1L [EA(L-)2 + 1EI ()2 + EAh(u (a-)
2 0 ax 3 ax ax ax

+ GA42 + 2cGAb(a ) + cGA()2 dx

T r1 r2h + h3- uh dx
T = 2pB u2h + w2h + 2 + uh2 dx
2 01 3 L




Substitution of Eqns. 3.10 into Eqns. 3.25, 3.26, and

Eqns. 3.25, 3.26 into Eqns. 3.20, 3.21 will result in a

stiffness matrix the same as Eq. 3.11 and a mass matrix as



pB x

3 i M dx

2h2I M Kj dx

+2 jKiMj dx

hfjK Kj dx


By substituting the shape functions in Eqns. 3.12

into Eq. 3.27 we then finally obtain the mass matrix as

pBL x


0 h










3.2.3 Lumped mass matrix

For some dynamic problems, use of lumped (or

diagonal) mass matrix is of considerable computational

convenience. To be more specific, if central difference

method is used in the dynamic finite element analysis the

calculation of the inverse of mass matrix can be avoided for

each time step. We now derive the procedure of the lumped

mass matrix as follows.

If numerical integration is used to evaluate the mass

matrix, we shall write a typical term as a summation

m = Jfi (x)pfj (x)dA = Wq[f, (Xq)pfj (xq)] (3.29)
fA q-1

where q refers to the sampling point at which the integrand

has to be found, p is the total number of the sampling

points for each element, fi is the general shape function,

Wq gives the appropriate weighting, and A is the domain of


If an integration scheme which uses only nodal points

for sampling is devised and which possesses the same order

of integration, such lumping will not affect the convergence

rate [19], and then all except one shape function fi are


mij = 0 (i j)

and the matrix becomes diagonal.

For the displacements 0 and u, which are the only

variables at two end nodes, the corresponding mass terms of

each element are simply split up equally into two parts and

each of them is lumped at one end node with a corresponding

weighting 0.5. Since the shape function for w is quadratic,

we shall therefore seek an integration formula which

integrates exactly a quadratic polynomial f(x) = a, + a2x +

a3x2 ,

f(x)dx = W1f(0) + W2f( ) + W3f(L)


where Wi represents the weighting of node i. Eq. 3.30 gives

W1 = W3 = L/6 and W2 = 2L/3. Using above procedure we obtain

the lumped mass matrix for beam element as

pBL x


3.3 Performance of Offset Beam Finite Element

In this section we examine the performance of the

offset beam element by solving some example problems, for

which the analytical solutions are available for comparison.

A Timoshenko cantilever aluminum beam is used with geometric

configuration shown in Fig. 3.2. The material constants and

dimensions are

E = 70 GPa G = 26.9 GPa L = 20 cm

B = 20 cm h = 1 cm P = 2.5 kN

a = 10 cm

where E and G are the Young's modulus and shear modulus, L,

B, and h are the length, the width, and the thickness of the

beam respectively, P is the magnitude of the applied load,

and a is the length of end notch.

3.3.1 Example 1: Displacement and rotation of tip of the

For a Timoshenko cantilever beam with a load P at the

end and with no crack we have

w = x + 2Px2 + -x (3.32)
6E1 2E1 icGA

2EIx Px (3.33)
2E1 El

and Table 3.1 shows the comparison between the finite

element solution and the analytical solution including the

relative error, as well as the convergence of the finite

element results.

Figure 3.2 Configuration of a sample cantilever beam

Table 3.1 The solution of the displacement and the rotation
at the tip of the beam by analytical method (Eqns. 3.32 and
3.33) and beam element method.

Number of Analytical
elements 4 8 16 solution

Beam tip
displacement(cm) 0.5368 0.5638 0.5710 0.5725

Error (%) 6.23 1.53 0.27 N/A

Beam tip
rotation -0.04286 -0.04287 -0.04291 -0.04286

Error (%) 0.00 0.0003 0.114 N/A

3.3.2 Example 2: Strain energy release rate of an end
notched beam

In the second example we consider an end notched

cantilever beam as a candidate to solve for the static

strain energy release rate of Mode II type fracture, GII, by

the offset beam element method and the analytical method.

The offset function of beam element is to be demonstrated.

As we will see in next chapter, the strain energy

release rate G can be expressed as

P2 ac
G a
2B aa


where C is the compliance of the beam. Since C is equal to

w,/P and wp is the displacement at the location where the

load is applied, then we have

G PWp (3.35)
2B 8a

For a cantilever beam with a long crack which runs in

the middle of the beam and for which the material property

is symmetric about the crack (Fig. 3.3a), the contact force

has been proved to be equal to half of the load and

distributed over a very small area right underneath the

location where the load is applied [20]. Therefore the beam

with a load on the top surface at the beam tip will

demonstrate a very similar deformation to the beam with two

loads, each of which is half of the original load and is

applied separately on the top or bottom half of the beam

(see Fig. 3.3b).

Since the crack is long, the top and the bottom parts

of the beam at the cracked area deform just like two beams

with the same load and configuration. Therefore the

deformation of the beam with long crack can be simulated by

the deformation of a beam with a varied cross-section, whose

bending rigidity changes from El of uncracked part to EI/4

of cracked part (see Fig. 3.3c). We now can easily derive

the tip displacement w, with the aid of classical beam

theory as

W = E( L3 + a3) (3.36)

Figure 3.3 Equivalent cases for calculating displacement and
strain energy release rate

Thus the closed form solution for G for a beam with a

long crack in the middle is available upon the substitution

of Eq. 3.36 into Eq. 3.35

3P2 a2
G 2E (3.37)

It may be noticed that the extension of crack da will

not change the cross-section area and hence the shear

deformation will not have any effect on G. Therefore, Euler-

Bernoulli beam theory, on which Eq. 3.36 is based, plays the

same role as Timoshenko beam theory.

For finite element method, Eq. 3.35 becomes

G P= A- (3.38)
2B Aa

Table 3.2 The solutions of the strain energy release rate
by analytical method (Eq. 3.37) and beam element method with
crack extension length fixed (Aa = 0.001 cm)

Number of Analytical
element 4 8 16 solution

G(kN/cm) 0.003013 0.003768 0.003963 0.004018

Error % 25.0 6.22 1.37 N/A

Table 3.3 The solutions of the strain energy release rate
by analytical method (Eq. 3.37) and beam element method with
number of elements fixed (16)

Crack extension Analytical
Aa (cm) 1 0.1 0.01 0.001 solution

G(kN/cm) 0.004335 0.003999 0.003967 0.003963 0.004018

Error % 7.89 0.463 1.27 1.37 N/A

Tables 3.2 and 3.3 illustrate the comparisons between

the analytical solution of Eq. 3.37 and beam element

solutions. In the meantime a two-dimensional finite element

analysis (Eight-Noded Isoparametric Two-dimensional

Quadrilateral Parabolic Element) is performed to solve for G

for the same problem. The G calculated by this finite

element method is 0.003750 kN/cm with Aa = 0.001 cm, 160

elements and the aid of Eq. 3.38. We then conclude that the

offset beam element gives convergent and considerably

accurate results for the problems of beams with interlaminar


3.3.3 Example 3: Natural frequencies of the beam

Since the dynamic crack propagation is the major

research motivation of this study, the dynamic response and

the natural frequencies of the beam are necessary for

examination. Unfortunately an analytical solution of dynamic

response for a Timoshenko cantilever beam is not available,

so at the present case only the natural frequencies are

investigated. In the next chapter, the dynamic crack tip

displacement as a function of time is calculated by the

offset beam element method and compared with that obtained

by two dimensional finite element analysis performed by

Grady and Sun [6].

The equation of the natural frequencies for an Euler-

Bernoulli beam is given as follows [21]:

.n = (P.L)2JEI/pL4

where EI, p, and L are the bending rigidity, density, and

length of the beam respectively, and (pL)2, determined by

beam characteristic equation, are given as 3.52, 22.4, 61.7

for the first three fundamental modes of a cantilever beam


It is very difficult to obtain the characteristic

equation for a cantilever beam by Timoshenko beam theory,

but rather easy for a simply supported beam. By comparing

the frequencies of a simply supported beam obtained from

Timoshenko beam theory with that from Euler-Bernoulli beam

theory, which will not be shown here, we found the

difference to be negligible. It has to be noticed that if

composites are used, whose G is comparatively small, the

natural frequencies of an Euler-Bernoulli beam can be very

different from those of Timoshenko beam.

Table 3.4 First three natural frequencies of vibration of
Euler-Bernoulli beam solution and beam element solution.

Number of Analytical
element 4 8 16 solution

Mode 0.1381 0.1330 0.1317 0.1318

Error % 4.78 0.910 0.0759 N/A

Mode 1.3693 0.8882 0.8312 0.8388

Error % 63.2 5.89 0.906 N/A

Mode 4.1806 2.7821 2.3512 2.3105

Error % 80.9 20.4 1.76 N/A

Table 3.4 shows the good agreement between finite

element results and analytical results as well as the high

convergence rate of the finite element solution.

Next, we compute the natural frequencies of top and

bottom beams joined by a distributed foundation spring as to

be explained and discussed in the next chapter. This

foundation spring model will be used in modeling the

interlaminar crack propagation. For finite element analysis,

we apply discrete springs to join the nodes of the pertinent

top and bottom beam elements. The constants of distributed

springs used here are EB/h for w, which is given by Kanninen

[22], and GB/h for u. Table 3.5 shows the agreement and

convergence for natural frequencies of the beam with


Table 3.5 First three natural frequencies of vibration of
Euler-Bernoulli beam solution and the spring model beam
element solution.

Number of Euler-Bernoulli
element 4 8 16 beam solution

Mode 0.1379 0.1328 0.1317 0.1318

Error % 4.63 0.759 0.0759 N/A

Mode 1.3642 0.8826 0.8327 0.8388

Error % 62.6 5.22 0.727 N/A

Mode 4.1806 2.7498 2.3458 2.3105

Error % 80.9 19.0 1.53 N/A

Thus we can conclude that the one-dimensional offset

beam element has a satisfactory performance for static,

dynamic, as well as the beam fracture problems.

3.4 Newmark Method

The offset beam element will be used to solve the

dynamic delamination propagation problems in the next

chapter as expected. The dynamic finite element equilibrium

equation mathematically represents a system of linear

differential equations of second order with damping term


MU + KU = R


where M and K are mass and stiffness matrices, R is the

external load vector, and U and U are the displacement and

acceleration vectors of the finite element assemblage.

In practice, we apply a finite element scheme of

static equilibrium at time t which includes the effect of

inertia forces over the solution domain, and we apply finite

difference method (or called direct integration method) for

the time domain.

Newmark method is one of the direct integration

methods and has been most widely used. It is found to be

simple, accurate, and is guaranteed to be stable [23]. One

can consider Newmark method as an accelerated direct

integration method,

Ut+At Ut
At = (1 6)Ut + SUt+At (3.40)

(Ut+At Ut)/At Ut
( At = (0.5 a)Ut + aUt+At (3.41)

and MA U+A + KA U = 0 (3.42)
t+At t+At t+At t+At

where 6 and a are the acceleration parameters to be

determined to obtain the integration accuracy and stability.

Here 6 = 0.5 and a = 0.25 are selected as suggested by Bathe


By rewriting Eqns 3.39 3.41 we obtain the final

explicit expressions,

Ut+At = (Kt+At + aoMt+At)- Mt+At (a0Ut + a2Ut + a3Ut)

Ut+At = a0 (Ut+At U ) a2 t a3Ut (3.44)

Ut+At = Ut + a6Ut + a7Ut+At (3.45)

where ao = 1/aAt2, a2 = l/aAt, a3 = 1/2a 1, a6 = At(l -

6), a7 = 6At [23].

Finally Ut+At, Ut+At and Ut+At can be determined only

by Ut, Ut and Ut, and mass and stiffness matrices through,

in turn, the application of Eqns. 3.43, 3.44, and 3.45.



4.1 Introduction

Fracture mechanics is an engineering discipline which

is primary applied to relate the maximum permissible applied

loads acting upon a structural component to the size and

location of a crack in the component. It can also be used to

predict the rate at which a crack will approach a critical

size in fatigue or by environmental influences, and can be

used to determine the time-dependent stress (or strain)

field during a rapid crack propagation in a solid. Current

damage tolerance assessment procedures are now available

that make effective use of these capabilities for materials

with a nonhomogeneous and anisotropic property.

We consider fracture as the formation of new surfaces

in a solid medium in a thermodynamically irreversible

manner. The study of the fracture process for a given

structure, a composite structure for example, requires the

simultaneous consideration of such widely diverse factors as

the macroscopic effects, which include the environmental and

loading conditions, stress waves, stress states particularly

around macroscopic imperfections where the crack is likely

to initiate or propagate, and the microscopic phenomena,

which include fiber breakage, debonding between fiber and

matrix, matrix dislocation and anything that takes place at

the locations where the crack nucleates and grows.

Therefore, due to the highly complex nature of the

phenomenon and the consequent lack of a full physical

understanding as well as the lack of sufficiently powerful

mathematical tools, at the present time there is no

consistent single theory dealing with all the relevant

aspects of fracture, even if a homogeneous continuum is


For the composite materials, the precise formation

and growth mechanisms are not always well defined. It is

generally assumed, however, that small flaws exist in the

material prior to loading. These flaws may be a result of

residual stress due to the curing process, external impact

damage, environmental degradation, or the fabrication

process. At some level of loading, these small flaws grow

and join together to form a single dominant interlaminar

crack. In the presence of the crack, the principles of

fracture mechanics can then be applied.

The linear elastic fracture mechanics will be applied

to composite materials for solving the designated problem

with the following assumptions:

1) The material is homogeneous, orthotropic, exhibits a

linear elastic behavior and undergoes brittle fracture.

2) There is only one dominant crack and the size of it is

sufficiently large compared to the characteristic dimensions

of the microstructure.

3) The crack always lies in between two layers and grows

along the interlaminar layer in a self-similar manner.

4.2 Discussion of Energy Release Rate

In the macroscopic approach to the problem, it

becomes necessary to devise a model for the actual

phenomenon and postulate a criterion for fracture. Among

these available criteria, most of them belong to either the

category of maximum quantity, like the stress intensity

factor K, the crack tip opening displacement COD, or the

category of energy balance, like the J-integral and the

strain energy release rate G. The energy balance criteria,

as we have known, have been widely accepted in different

applications due to their generality, physical soundness and

flexibility. If energy as a whole is concerned, the analysis

of the complex stress states and microscopic phenomena at

the crack tip then can be avoided. Consequently, the

experimental measurement or the concentration of elements in

finite element calculation at the crack tip can also be


In essence, all the energy balance criteria are based

upon a simple thermodynamic notion that the crack will ensue

or continue to propagate if, for a unit increase in crack

surface, the increase of externally added or internally

released energy is greater than or equal to the amount of

stored, kinetic and dissipated energies. In other words, the

rates of input, stored, dissipated, and kinetic energies are

necessarily to be balanced.

The energy balance for a plane problem can be

expressed as

S U + T W + = 0 (4.1)
B da I

where U, T, W, and D are stored energy (elastic strain

energy), kinetic energy, input energy (potential energy of

the external load), and dissipated energy (mainly surface

energy for brittle materials) respectively, B is the

thickness of the material, and a is the crack length.

Further, Eq. 4.1 can be written as

ld 1 dD
1 d (W U T) B da (4.2)
B da B da

We then define the expression in left hand side as

energy release rate G (or called crack extension force), and

the right hand side as the critical energy release rate Gc

(or the crack resistance). Gc is a constant for the given

material initially, but may depend on the crack velocity.

When G < Gc the system keeps an energy balance and no crack

extension happens. G = GC is the critical point for crack

propagation. If it takes place, it turns to a new energy


The simplest mathematical formula for G in static

case is conceived as

1 d
G d(W U)
B da
1 d 1
B da (q 2Pq)

SP 2 PB (4.3)
2B da

where P is the external load and q, is the generalized

displacement at the location where the load is applied.

To be more general, one can rewrite Eq. 4.2 over a

discrete solution domain as

G da(W U T)
B da

1 d
VB dt(W U T)

1id 1 1
VB -dtj(Pq X jqIj q-qj)
B ij ij

-v(PqP j jqqj CXXl j qj qj
VB ij ij

-1 i j XXM j q ~)
ij ij

where V is the crack extension velocity, Kij and Mij are the

components of stiffness and mass matrices respectively, and

the quantities with dots at the top are their derivatives

with respect to time.

Consider the following expression, which is the part

of the expression above,

Pqp XXKI j ,q XMj j 4 q"
ij ij

= XPqI 6P jCI Kj j qj XMj j q4 q4
i ij ij

= Xq (P6~ XK,,q, M~,j j)
i j s

The expression within the parenthesis is identically

zero by the definition of dynamic equilibrium. Therefore the

expression for G becomes

G = -B X(Ki j qq + Mi qjg) (4.4)

Equation 4.4 derived is comprehensive and fundamental

in computing G by finite element method. However, the crack

growth velocity V can not be obtained a priori, and

technically taking derivatives of the stiffness and mass

matrices is not easy. In the next section we will discuss

several other methods which are more appropriate to the

present study of the interlaminar crack propagation.

4.3 J-Integral and Crack Closure Integral

Dynamic fracture mechanics is for those fracture

problems in which not only inertia forces are included but

the loading and the crack size vary dynamically. The

involvement of system kinetic energy, crack growth velocity,

and the contact of the crack surfaces (which causes

friction) due to impact, etc. make the problem rather

complex. For a beam laminate with an embedded delamination,

the zero-volume integral proposed by Farris and Doyle [24]

can be applied. This integral is a special case of J-

integral, which is equivalent to the energy release rate in

linear elastic fracture mechanics. It expresses G as a

function of the bending moments, shear and axial force

resultants of top and bottom parts of the beam at the crack

tip cross-section.

Since composites are very brittle materials, the

plastic zone at the crack tip is negligible. The path of the

J-integral finally shrinks to one line perpendicular to the

crack at its tip. So it is also called zero-volume integral.

If the crack is in the middle of the beam, we have

h3 [ M ]2 (.
S- 72E I I (4.5)

where M, I, E, and h are the bending moment, the moment of

inertia, the Young's modulus, and the thickness of the beam

respectively. The quantities with bars at the top stand for

those of sublaminates.

Irwin's crack closure integral [25] has been

successfully applied to the beams with embedded delamination

[6]. The theory simply states that if the crack extends by

Aa, the energy absorbed in the process is equal to the work

required to close the crack to its original geometry. Using

a polar coordinate system with the origin at the extended

crack tip, this statement in equation form is as follows:

1 Aa
G = lim 2AaJ {y (Aa-r,0)Av(r,7r) + Tr, (Aa-r,0)Au(r,7)}dr

where Av and Au are the relative sliding and opening

displacements between corresponding points on the crack


In finite element modeling, the stresses over the

crack extension Aa are replaced by the nodal forces of a

whole element which is of the same size as the unit crack

extension. Then Irwin's theory states that energy release

rate is equal to the work required to close two nodes at the

crack tip. If Aa is small enough and xi+I- xi = xi- x.i1,

the forces required to close the crack can be replaced by

the forces of the crack tip node i, and the displacements at

the node i+l can be considered as the displacements through

which the forces do the work [26] (see Fig. 4.la). For beam

element formulation, the energy release rate G is

G 2Ba[aQiAwi+1 + MiAi+1 + NiAui+1] (4.6)

where Q, M, and N stand for the forces corresponding to the

displacements w, 0 and u, and Aw, AO and Au are the relative

displacements of upper and lower crack surfaces.

For the work done in transverse direction, the part

of the middle node must be included. i.e., G is the sum of

the terms in the bracket of Eq. 4.6 and the last term

QjAWj+,, where j is the middle node between nodes i and i-1,

and j+l is the middle node between nodes i and i+l (see Fig.

4.la). Finite element analysis has shown that Qj is always

negligibly small and so is the last term.

In order to compare the theoretical modeling with

experiment results, we select the specimen used in Ref. 6 as

the candidate, which is 20-ply [90/0]ss T-300/934

graphite/epoxy laminates. Its configuration is in Fig. 4.2

and the effective material constants are as follows:

EI = 72.4 GPa E2 = 10.3 GPa

v12 = 0.33 p = 1.58x10-5N-s2/cm4

The impact load history is also given by Grady and

Sun [6] as

S S i~I

Rigid Element

or Spring Element

Figure 4.1 Finite element nodes at crack tip

i+1 j+

i L

- -


S155 m/s


0.173 J

Unit: m

Figure 4.2 Dimensions of the specimen used in Ref. 6

F F sin(nt/T) t < T
F(t) =
S0 t T

where the amplitude of force history Fo and the contact

duration T are 1530 N and 125 ps.

Based upon the data given by Grady and Sun [6], the

energy release rate G is computed by the offset beam

element. The number of elements in this case is 1000 and the

time step is lps. Figures 4.3 and 4.4 show the results of

energy release rate G as a function of time obtained by J-

integral and crack closure integral respectively. From Fig.

4.8 it can seen that these two methods agree well. The

relative error does not exceed 2.8% at the time period

concerned. Figure 4.5 also shows the agreement of the

deflections of the crack tip node as a function of time

computed using the present element and two dimensional

element in Ref. 6.

4.4 Rigid Element

As mentioned before, the element used in this study

for finite element analysis consists of top and bottom

elements and each pair of them is connected by three nodes.

In order to simulate the crack growth the node at the crack

tip will break into two, which causes some computational

difficulties in remodeling the finite element mesh. But this

120 1
Time (micro-second)

Figure 4.3 Energy release rate by zero volume integral

0 40 80 120 160 200 240
Time (micro-second)

Figure 4.4 Energy release rate by crack closure integral

0.0017 -
0.0016 -
0.0015 A Ref. 6
0.0014- this study
0.0013 -
0.0012 -
0.0011 -
rE 0.001 -

C 0.0006 -
S 0.000 displacement
S 0.0007 -
0.0006 -
0.0005 -
0.0004 -
0.0003 -
0.0002 -
0.0001 -
0 -A
-0.0001 -
0 40 80 120 160 200 240

Time (micro-second)

Figure 4.5 Transverse displacement at left crack tip

problem can be solved by applying the so called "rigid"

element, which ensures the continuity of transverse

displacement, axial displacement, and rotation of top and

bottom elements. In other words, top and bottom elements

connected by a dimensionless rigid element will function in

the same way as if they were joined by identical nodes.

Removing the rigid element is equivalent to creating crack

or crack growth (Fig. 4.1b).

Each rigid element has nine degrees of freedom. In

addition to the three displacements each at both top and

bottom nodes, the three generalized forces that the rigid

element transmits are also considered as unknown degrees of

freedom. The constraint equations corresponding to the rigid

element can be written as follows:

F 0 0 -I q

F2 = 0 0 I q (4.7)

0 -I I 0 F

where {F}1 and {F2), {q ) and {q2}, are the nodal force

vectors and displacement vectors respectively, at node 1 and

2, corresponding to the top and the bottom elements. {F} is

the vector of forces transmitted by the rigid element and is

also considered as an unknown. [I] is the 3 x 3 identity

matrix. This formulation can be understood as the Lagrange

Multiplier used to satisfy the equality constraints.

Another advantage for the rigid element is that the

forces between two corresponding nodes can be obtained

directly. It is especially convenient if the contact force

is needed for determining the friction force.

The coefficients in Eq. 4.7 should be multiplied by a

suitable factor, say Kf, to be of the same order of the

magnitude as that of the element of the stiffness matrix, in

order to avoid introducing a large round-off error. Thus Eq.

4.7 becomes

F 1 0 0 -Kf q

F2 = 0 Kf I q (4.8)

0 -KI KfI 0 F/K,

In the delaminated portion we assume that there is no

friction between the contact surfaces. Thus we need the

rigid element that transmits only compressive force in the

transverse direction. Then the stiffness matrix will be of

size 3 x 3 as shown below,

{ I' 0 0 -Kf wI

Fw = 0 0 Kf Wz (4.9)

0 -Kf Kf 0 F,/Kf

where F, is the component of the vector {F} corresponding to

the transverse displacement w. It may be noted that in the

delaminated portion, the rigid element will be effective

only when the crack surfaces are in contact or overlapped,

i.e. wl w2 > 0. If wi w2 : 0, then the rigid element

will not be connected and consequently a gap will be formed

as it was supposed to be.

For the middle node of the beam element, which has

only one degree of freedom in transverse direction, the

constraint equations of the rigid element will be the same

as the one used in the delaminated portion described in Eq.


It has been noticed that the global stiffness matrix

will not be diagonally dominant if the matrices of the rigid

elements are assembled into it. Especially, some linear

system solvers can not solve the system with zero magnitude

elements on the diagonal. Therefore the application of an

IMSL subroutine as a linear system solver is suggested since

it has the facility to switch rows and columns to make the

matrix diagonally dominant and then to solve it.

4.5 Elastic Foundation Spring Model

The so called spring model was first proposed by

Prandtl [27,28], which was based on the principle of

breakage of nonlinear atomic bond at crack tip. The elastic

foundation spring model, developed by Kanninen, has a

distinctly different physical interpretation. His spring is

ideal, linearly elastic, and continuously distributed.

4.5.1 Principle of energy balance

Kanninen and his co-workers [29,30,31] have developed

various elastic foundation spring models for dynamic

analysis of unstable crack propagation and arrest in Double

Cantilever Beam (DCB) specimens. The fundamental advantage

of the method lies in the simplicity in calculating the

dynamic energy release rate G, especially with the

application of finite element analysis. However, the

difficulty is to define the proper spring constants for a

given delaminated beam. In this part the attention is

focused on the application of the spring model to the finite

element analysis for the problem concerned.

The spring model used here extends the application of

the energy balance principle. Equation 4.1 expresses the

balance of elastic strain energy U, kinetic energy T, work

done by external load W, and dissipated energy D for a given

system. For brittle materials like composites, dissipated

energy is almost equal to crack surface energy since the

plastic work and viscous dissipation are negligible. For the

spring model, the crack extension is simulated by the

rupture of springs, which is equivalent to the formation of

new crack surface. Thus if a crack extension da is assumed,

for a beam with spring the total spring energy U* will

decrease by dU*, whereas in the actual beam, the total

surface energy D will increase by dD. In principle, the

distributions of stress, strain and displacement, and

furthermore the potential energy, the kinetic energy, and

the work done by external loads, will remain the same if a

beam with crack is replaced by two separate beams partially

(at the uncracked part) joined by a distributed spring. Thus

we can conclude that dU* is equal to dW and Eq. 4.1 then


Sd U + T W + U = 0

and consequently we obtain

1 d 1 dU*
G da(W U T) Bda (4.10)

Hence for the designated problem the spring model takes the


G = 2B [C (Aw)2 + C2 (A)2 + C3 )2] (4.11)

where C1, C2, and C3 are the spring constants for w, 0, and

u respectively.

It also has to be noticed that dU* = dW only when the

crack extension occurs. In fact, (1/B)dU*/da is equal to G

while (1/B)dW/da is equal to Go, the critical value of G.

4.5.2 Investigation of spring constants

In Ref. 29 Eq. 4.10 is derived in a different way. By

using Timoshenko beam theory and assuming lower half of the

top part of the DCB specimen (because of symmetry for Mode I

cracking, only the top part is used) is elastic and the

upper half is rigid, Kanninen derived a formula for the

energy release rate for Mode I fracture as

G =[ 2E w2 (4.12)

where E, B, and h are the Young's modulus, the width, and

the thickness of the beam, and w is the transverse

displacement at the crack tip.

Equation 4.12 yields that the spring constant for w

is 2EB/h. The spring constant will be the half, i.e. EB/h,

if a whole DCB specimen is considered. Apparently, this

constant comes from the assumption that half of the top part

of the beam behaves in an elastic manner, and the other half

is rigid.

By definition, beams are entirely rigid in the

transverse direction. Returning to the original model, i.e.

a beam with a distributed elastic foundation spring on the

bottom which can be understood as half of a double

cantilever beam (DCB), one may observe the equilibrium

equation as

EI + kw = 0 (4.13)

where El is the bending rigidity of the beam, k is the

spring constant, and w is the transverse displacement. For

simplicity an Euler-Bernoulli beam is used. The loading

condition is set to produce the pure Mode I type fracture

and only transverse spring corresponding to w is included.

In order to compare, this condition is extended to the

finite element analysis. The material constants used are the

ones introduced in section 3.3. With the appropriate

boundary conditions which lead to a pure Mode I type

fracture, Eq. 4.13 is solved for w. By equating the energy

release rate of the spring model, which is a function of

spring constant k, to the one obtained by J-integral, one

then determines the appropriate spring constant.

Surprisingly, the constant is found to be infinity. In other

words, the energy release rate of the spring model goes to

that of J-integral asymptotically as the spring constant

goes to infinity, as shown in Fig. 4.6. The detailed

derivation appears in appendix B. This result would seem

less surprising if one understands that any non-infinite

spring constant will result in a non-zero crack tip

displacement, which does not precisely describe the physical



Spring Model (Discrete Spring)
0.7 Spring Interval o 0.625 cm
J integral A 0.2 cm
x 0.04 cm
0.6 -


0 0.4

\ Spring Model
S0.3 1(Continuous Spring)
0 .3



0 O-x

-3 -1 1 3 5 7


K: The spring constant for trial

K1: The spring constant used in Ref. [29]

Figure 4.6 Comparison of the results for energy release rates by J-integral,

distributed spring model, and discrete spring model for Mode I fracture

phenomenon. As a matter of fact, the spring constant

suggested by Kanninen, though not perfect, is a good enough

approximation (Fig. 4.6). It is worth pointing out that if a

given spring constant is large enough, its variation may

produce a corresponding variation of the displacement w at

the crack tip but the energy release rate G will not greatly


Unlike the distributed spring discussed above, the

discrete spring, used commonly in finite element analysis,

has a different physical interpretation. Figure 4.6 shows

that as the spring constant increases, the energy release

rates of finite element calculation do not asymptotically go

to the J-integral result but to zero. This can be explained

as follows. If the spring constant is considerably large,

the spring model is comparable with the crack closure

integral. The crack closure integral discussed in section

4.3 states that for finite element analysis the energy

release rate is equal to the work required to close two

nodes at the crack tip (see Eq. 4.6). If the displacement

and the force are understood as the spring displacement and

the spring force in Eq. 4.6, the crack closure integral

model can be interpreted as the spring model. Nevertheless,

the displacement at node i+l is always larger than the

displacement at node i, unless Aa goes to zero. This

explains why the energy release rate goes to zero for

discrete spring, instead of to the right answer, if spring

constant increases. Therefore, in order to apply the spring

model to finite element analysis, one needs to loosen the

discrete springs from infinitely stiff by multiplying a

coefficient ( to the Kanninen's spring constant (here

consider Mode I type fracture and transverse displacement w

only). This coefficient is selected by comparing the energy

release rate obtained by spring model with the one by J-

integral. Table 4.1 shows that ( becomes larger as the size

of element is reduced. It is also shown in Fig. 4.6 that

intersection of the G- curve and the J-integral result

moves toward right as the size of the element is reduced. If

the size of element approaches zero, the coefficient of

spring constant ( will asymptotically approach to infinity.

In Eq. 4.6 one finds that Aw/Aa will remain a

constant irrespective of the change in size of element Aa,

since the force at the nodal points will not change greatly.

By comparing Eq. 4.6 with Eq. 4.11, one may realize that for

spring model the coefficient ( multiplied by (Aa)2 will

become a constant not greatly varying with Aa. In other

words, the energy release rate can be expressed as

G = j[j(Aw)2 = C )
2 h 2h Aa

where EB/h is the Kanninen's spring constant for a whole DCB

specimen. Apparently, as Aa becomes small enough, C = e(Aa)2

is independent from Aa and is dependent only on material

constants and geometric configuration. Table 4.1 supported

this concept and it can be used to determine the coefficient

( for spring model.

Table 4.1 Coefficient ( of Kanninen's spring
the different size of element

constant for

Both static and dynamic examples have proved that if

an appropriate spring constant is selected, the discrete

spring model gives an identical result to those from J-

integral and crack closure integral. For the same load,

material constants, and dimensions introduced in section

4.3, the dynamic energy release rate is also calculated by

spring model (Fig. 4.7) with the choice of the coeeficient

of spring constant as 9.0. Figure 4.8 illustrates that the

No. of Element Size Coefficient ( C=(Aa)2
Elements (cm)

100 0.40000 0.9506 0.15210

200 0.20000 4.2658 0.17063

300 0.13333 10.000 0.17778

400 0.10000 18.197 0.18197

500 0.08000 28.840 0.18458

600 0.06667 41.976 0.18156

700 0.05714 57.544 0.18788

800 0.05000 75.683 0.18921

900 0.04444 96.161 0.18994

1000 0.04000 119.124 0.19060

0 40 80 120 160 200 240
Time (micro-second)

Figure 4.7 Energy release rate by spring model

0 40 80 120 160 200 240
Time (micro-second)

Figure 4.8 Energy release rates by three different approaches

agreement of the results from three methods is satisfactory

based upon the same number of elements and time step.

4.5.3 Dynamic crack propagation

Whenever the energy release rate G reaches its

critical value G,, the crack starts to propagate. If the

spring model is applied to finite element analysis, the

delaminated composite beam performs simply as two parallel

and separate beams joined together with a number of springs

at the uncracked portion and these springs break one by one

to simulate the crack propagation. Apparently, the geometric

configuration and the position and the size of the crack of

the specimen, as well as the shape and the velocity of

impactor can significantly affect crack propagation. In

fact, if the tip of the embedded crack is not sufficiently

distant from the impact site, the influence of the

indentation and the distribution of the contact stress must

be considered. But for the designated problem, the error

introduced by neglecting the contact behavior should be

small [6].

Based upon the same critical energy release rate Gc

(350 N/m) and the loading history F(t) given by Grady and

Sun [6], the crack propagation in a delaminated composite

beam subjected to impact is simulated by the spring model in

finite element analysis. Actually if the crack growth

occurs, the critical energy release rate Gc is a function of

the crack growth velocity (a linear function usually, if the

velocity is not too high), which is larger than the critical

energy release rate without crack growth [29]. The finite

element computation for crack propagation has shown that the

involvement of the velocity coefficient for Gc does not

significantly change the curve of crack tip position versus

time. Therefore it is not considered in the current study.

Figure 4.9 illustrates an acceptable agreement

between the finite element results of the present study and

the experimental results from Ref. 6 for the crack tip

position as a function of time. In computing the dynamic

energy release rate G, At is assigned as 1 As, Aa is

1.37x10-3m which is the length of element, the total number

of elements is 256, and the coefficient selected for spring

constant is 1.8. G is calculated at each time step. If G is

smaller than G,, the crack does not propagate. If G is equal

to or greater than Gc, the crack starts to propagate with a

increment Aa. However, if G is much larger than G,, it

implies the maximum crack propagating velocity assumed,

Aa/At, is not appropriate. In essence, if a larger maximum

velocity is assumed, after a crack extension Aa it may take

a few more (e.g. n) time steps to produce a new crack

extension. The current crack growth velocity becomes

Aa/(nAt). Nevertheless, if a smaller maximum velocity is

assumed and G is much larger than Gc, a larger G, is wrongly

assumed and the results will be inaccurate. So we would


0.024 -
158 m/s
0.022 a

0.02 -

0.018 -

0.016 -

C Experimental Data'from Ref. 6
A 0.012 -
1 0.01 --- FEM Results by This Study

0 0.008

0.006 -

0.004 -

0.002 -

0 --- m -- -- -l--mm--------m---
0 200 400 600 800
Time (micro-second)

Figure 4.9 Experimental and finite element results of left crack
tip position: crack extension versus time

rather assume a larger maximum velocity, either to reduce At

or to increase Aa for adjustment. Usually for the

convenience of programming, both Aa and At are fixed. We

need to estimate the maximum crack growth velocity and

consequently to choose the corresponding Aa and At.

4.6 Buckling Effect

Delamination is a cause of progressive stiffness loss

and subsequent strength degradation. In the case of an

unsymmetric delamination layer, the postbuckling deformation

and the strain energy release rate are strongly affected by

bending-compressing [32] and shear-compressing coupling. It

has been shown that the damage tolerance is determined

either by the material fracture toughness or the buckling

load of the delaminated sublaminate, depending on the size

and location of the interlaminar crack. Larger

delaminations, and those closer to the outer surface of the

composite beam were shown to be more susceptible to

delamination buckling [33]. The calculation of Mode I and

Mode II strain energy release rates showed that for a beam

with a large delamination and under the transverse impact

loading, although GII reaches significantly greater

magnitude than GI, the magnitude of GI controls the cyclic

growth rates of the delaminations [34]. In fact if the

delamination is large enough it causes buckling, followed

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs