DYNAMIC DELAMINATION PROPAGATION
IN COMPOSITE BEAMS UNDER IMPACT
By
SHOUFENG HU
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1990
ACKNOWLEDGEMENTS
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
possible.
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.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS ...................................... ii
LIST OF FIGURES ...................................... v
ABSTRACT .......... ................................ vii
CHAPTERS
1 INTRODUCTION ...................... .......... 1
1.1 Background ................................... 2
1.2 Objective of This Study ....................... 4
2 ANALYSIS OF INTERLAMINAR SHEAR STRESS FOR AN
ORTHOTROPIC BEAM DUE TO STATIC INDENTATION .... 6
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 Loadcontact length relations ........... 22
2.4.3 Shear stress investigation ............... 23
2.4 Conclusions ................................... 29
3 OFFSET BEAM FINITE ELEMENT FOR STATIC AND
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
iii
3.4 Newmark Method ................................ 58
4 ON FRACTURE BEHAVIOR OF A DELAMINATED
COMPOSITE BEAM SUBJECTED TO IMPACT ........... 61
4.1 Introdution ................................... 61
4.2 Discussion of Energy Release Rate ............ 63
4.3 JIntegral 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 CONCLUSIONS AND FUTURE RESEARCH EXPECTATIONS .. 113
6.1 Conclusions ............. .................. .. 113
6.2 Future Research Expectations .................. 115
APPENDICES
A DERIVATION OF STRESSES BY SUPERPOSITION
METHOD .......... ............................. 118
B DETERMINATION OF SPRING CONSTANT FOR
SPRING MODEL ...... ........................... 120
C ON THE BUCKLING EFFECT OF CRACKED LAMINATES ... 125
D THEORETICAL PREDICTION OF IMPACT FORCE
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
LIST OF FIGURES
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) halfplane solution,
(b) Timoshenko elastic beam solution, and
(c) quasielastic 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
DYNAMIC DELAMINATION PROPAGATION
IN COMPOSITE BEAMS UNDER IMPACT
By
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 nodeshifting technique is used to solve two
dimensional fracture problems using a onedimensional beam
element. The technique has been successful in calculating
the energy release rate of a beam with an embedded
vii
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 Jintegral, 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.
viii
CHAPTER 1
INTRODUCTION
Improved specific strength and stiffness properties
have made fiberreinforced 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 fibermatrix 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
5
theoretical models. Chapter 6 includes the summary of the
present study, concluding remarks, and suggestions for
future studies.
CHAPTER 2
ANALYSIS OF INTERLAMINAR SHEAR STRESS
IN AN ORTHOTROPIC BEAM DUE TO STATIC INDENTATION
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
distribution.
A simple, yet accurate procedure is developed in the
present study for obtaining an analytical solution for
stresses in an orthotropic beam under threepoint 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 halfplane 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
observed.
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 halfplane, an elastic beam suggested by
Timoshenko and a quasielasticity 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 halfplane 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
y
(a)
0.1
6.25
1: h
 Cross Sections
for Normal Stresses
.......... Cross Sections
for Shear stress
Unit: mm
(b)
Figure 2.1 A homogeneous and orthotropic beam subject to a uniformly
distributed load over a length 2a
wI
*C~l~
I
v
I
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 halfplane can be calculated
from the point force solution as
a (x,y) = (P/2a) gi (xC,y)dC, i = 1,2,3 (2.3)
a
where the superscript h stands for halfplane.
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 halfplane 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 quasielasticity
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
(a)
I 1
(b) (c)
Figure 2.2 Superposition of (a) halfplane solution, (b) Timoshenko elastic
beam solution, (c) quasielastic solution
1T[
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)
Lx
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 selfequilibrating 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 quasielasticity 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
equilibrium.
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 threepoint
bending then can be obtained by superposing stress again as
5
+0.5
Figure 2.3 Superposition of the three cases for stresses of a short beam
due to three point bending with distributed load and supports
Ill
 ...
ai (x,y) = Lo (x,y) + 0.5aL (xd, hy)
+ 0.5a4 (dx, hy) (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 nondimensionalized as
follows, in order to be compared with the stresses obtained
by whitney [7]:
Bh2 2Ba Bh
axx 2Pdxx yy p ayy xy prxy
2.3 Numerical Results
Whitney [7] performed a twodimensional Fourier
series elasticity analysis for an orthotropic beam under
threepoint (and fourpoint) 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 D2344 [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 (EightNoded Isoparametric
Twodimensional 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 closedform
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
x FEM
1.8
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
0 2 4
Yaxis (mm)
Figure 2.4 Normalized transverse shear stress Tau(xy) distributions at
x = 0.05, 6.25, and 12.0 mm
2
0
1 
2 2 
0X0 11
E 3 xmm x=12.5mm
D 4
c 5 +
6 5
Ref.7
7
.A FEM
cc 8
z 9 This Study
10 
11 
12
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
0.8
S FEM
0.7
 This Study
D 0.8
~0 o.
0.4
Z x=O. mm x=12.5 mm
Z 0.2
0 II
0.1
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
study.
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 halfplane.
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 halfplane
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 Loadcontact length relations
Similar to the situation discussed above, for a
relatively small radius indentor (or a thick beam), the
loadcontact length relation can be derived to be the same
as that for an orthotropic halfplane [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
[13]
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
composites
Type Material El E2 V12 G12
GPa GPa GPa
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
8
7
Graphite/Epoxy
6 + Glass/Epoxy
S5 Boron/Epoxy
+
E 4
3 
+ +
2
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.5
1.4
1.3 0 Graphite/Epoxy
1.2 
+ Glass/Epoxy
S 11 
1 Boron/Epoxy
S0.9 
0
C 0.8 
E
E 0.7 
M 0.6 
S 0.5 
0.4
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
investigation.
9
8
0 Uniform Loading
A Elliptical Loading
e Beam Material: Graphite/Epoxy
Load: 50 N
S 5 Contact Length: 0.4206 mm
(D 4
3
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
stress.
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.
CHAPTER 3
OFFSET BEAM FINITE ELEMENT
FOR STATIC AND DYNAMIC PROBLEMS
3.1 Introduction
A onedimensional 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 nodeshifted 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 xyz 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 twodimensional 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 twodimensional 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 yaxis, and w(x) as the displacement of the
beam in zdirection (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 crosssectional plane is the same,
i.e., there is no thickness stretch. These are two
fundamental assumptions for beam theory.
Z
X wu w wIu u
T T
Figure 3.1 Configuration of top and bottom beam elements with
nodes shifted
Using the straindisplacement and stressstrain
relations, we readily obtain
S= Ez dx (3.3)
dxdw
rxz = IG(O + d) (3.4)
where E is the longitudinal Young's modulus and G is the
shear modulus in the xz 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
threedimensional theory [16]. Cowper [17] provided a
detailed discussion of the determination of K. For a beam
with a rectangular crosssection 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
obtain
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 crosssection
area.
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 zdirection, is
dx + p(x) = 0
and it can be written in terms of displacements as
di GA( + )] + p = 0 (3.7)
where p(x) is the external distributed load in transverse
direction.
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 xdirection
dN
dx + t = 0
dx
where N is the normal force in longitudinal axis of the beam
and
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
(3.8)
and
dx EAd(u )] + t = 0 (3.9)
3.2.1 Method of weighted residualsa differential
formulation
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
obtain
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 ,
(3.10)
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 dMdx
fLdx f dx dx 2 dxd
fpNi dx
L
ftKidx
JL
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
(3.11)
From Eqns. 3.73.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
x
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
EA 0 EAh EA
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
EA
L
(3.13)
the corresponding displacements, in turn, are w 1 ul,
W2 W3, 3 3 and u3
3.2.2 Principle of total potential energya variational
formulation
A continuum problem usually has different but
equivalent formulationsa 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
yields
a8
0 i = 1,2,3,*,n (3.16)
Upon the substitution of Eq. 3.14 into Eq. 3.16 we
obtain
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)
ii
T = 2 i (3.19)
ii
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
(3.20)
(3.21)
Expressions for U and T in terms of nodal
displacements can be determined as
:Lh / 2 1
U = f (EIe2 + nGAy )dZdx
S Jh / 2
T i pB(u2 + u2)dzdx
0 J h/2
(3.22)
(3.23)
With the aid of straindisplacement 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
(3.24)
(3.25)
(3.26)
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
follows:
ShfNNjdx
JL
pB x
3 i M dx
2h2I M Kj dx
+2 jKiMj dx
hfjK Kj dx
JL
(3.27)
By substituting the shape functions in Eqns. 3.12
into Eq. 3.27 we then finally obtain the mass matrix as
pBL x
2h
15
0 h
15
h
30
0
0
h
15
12
12
h2
+
12
symmetric
h2
6
(3.28)
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 q1
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
element.
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
zero,
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)
(3.30)
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.31)
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
beam
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
(3.34)
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 crosssection, 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)
2BEI
It may be noticed that the extension of crack da will
not change the crosssection 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 twodimensional finite element
analysis (EightNoded Isoparametric Twodimensional
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
cracks.
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
[21].
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 EulerBernoulli 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 EulerBernoulli beam can be very
different from those of Timoshenko beam.
Table 3.4 First three natural frequencies of vibration of
EulerBernoulli beam solution and beam element solution.
Number of Analytical
element 4 8 16 solution
First
Mode 0.1381 0.1330 0.1317 0.1318
Error % 4.78 0.910 0.0759 N/A
Second
Mode 1.3693 0.8882 0.8312 0.8388
Error % 63.2 5.89 0.906 N/A
Third
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
springs.
Table 3.5 First three natural frequencies of vibration of
EulerBernoulli beam solution and the spring model beam
element solution.
Number of EulerBernoulli
element 4 8 16 beam solution
First
Mode 0.1379 0.1328 0.1317 0.1318
Error % 4.63 0.759 0.0759 N/A
Second
Mode 1.3642 0.8826 0.8327 0.8388
Error % 62.6 5.22 0.727 N/A
Third
Mode 4.1806 2.7498 2.3458 2.3105
Error % 80.9 19.0 1.53 N/A
Thus we can conclude that the onedimensional 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
neglected
MU + KU = R
(3.39)
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
[23].
By rewriting Eqns 3.39 3.41 we obtain the final
explicit expressions,
Ut+At = (Kt+At + aoMt+At) Mt+At (a0Ut + a2Ut + a3Ut)
(3.43)
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.
CHAPTER 4
ON FRACTURE BEHAVIOR
OF A DELAMINATED COMPOSITE BEAM SUBJECTED TO IMPACT
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 timedependent 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
assumed.
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 selfsimilar 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 Jintegral 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
averted.
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
balance.
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 qqj)
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 JIntegral 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 zerovolume 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 crosssection.
Since composites are very brittle materials, the
plastic zone at the crack tip is negligible. The path of the
Jintegral finally shrinks to one line perpendicular to the
crack at its tip. So it is also called zerovolume 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 (Aar,0)Av(r,7r) + Tr, (Aar,0)Au(r,7)}dr
where Av and Au are the relative sliding and opening
displacements between corresponding points on the crack
surfaces.
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
1
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 i1,
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 20ply [90/0]ss T300/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.58x105Ns2/cm4
The impact load history is also given by Grady and
Sun [6] as
S S i~I
0
Rigid Element
or Spring Element
Figure 4.1 Finite element nodes at crack tip
i+1 j+
i L
 
h
S155 m/s
0.129
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 (microsecond)
Figure 4.3 Energy release rate by zero volume integral
0 40 80 120 160 200 240
Time (microsecond)
Figure 4.4 Energy release rate by crack closure integral
0.0019
0.0018
0.0017 
0.0016 
0.0015 A Ref. 6
0.0014 this study
0.0013 
0.0012 
Impact
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 (microsecond)
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 roundoff 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.
4.9.
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 coworkers [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
becomes
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
form
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
d4w
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 EulerBernoulli 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 Jintegral, 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 Jintegral 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 noninfinite
spring constant will result in a nonzero crack tip
displacement, which does not precisely describe the physical
83
0.8
Spring Model (Discrete Spring)
0.7 Spring Interval o 0.625 cm
J integral A 0.2 cm
x 0.04 cm
0.6 
Ex
S0.5
0
0 0.4
\ Spring Model
S0.3 1(Continuous Spring)
0 .3
x
0.2
0.1
0 Ox
3 1 1 3 5 7
Log(K/K1)
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 Jintegral,
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
vary.
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 Jintegral 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 Jintegral 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 (microsecond)
Figure 4.7 Energy release rate by spring model
0 40 80 120 160 200 240
Time (microsecond)
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.37x103m 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.026
0.024 
158 m/s
0.022 a
0.02 
0.018 
0.016 
S0.014
C Experimental Data'from Ref. 6
A 0.012 
0
a
0
1 0.01  FEM Results by This Study
0 0.008
0.006 
0.004 
0.002 
0
0  m   lmmm
0 200 400 600 800
Time (microsecond)
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
bendingcompressing [32] and shearcompressing 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
