"2:53.: . x

 

 

 

I LIBRARY
Michigan State

I) 91):? 1/ Q W University

\

é

This is to certify that the
thesis entitled

MODELING AND DESIGN OPTIMIZATION OF GOLF
SHAFTS

presented by

NICHOLAS FRANK ABBRUZZESE

has been accepted towards fulfillment
of the requirements for the

MS. degree in Mechanical Engineering

 

 

Mam

Major Professor’s Signature
7/ 2 5'[c7 5‘

Date

MSU is an Affirmative Action/Equal Opportunity Institution

 

PLACE IN REIURN BOX to remove this chedtout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

NOV 0 3 ZOC3
3

 

vV 5‘ '03

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 mm.mts

 

MODELING AND DESIGN OPTIMIZATION OF

GOLF SHAFTS
By

Nicholas Frank Abbruzzese

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Mechanical Engineering

2005

ABSTRACT

MODELING AND DESIGN OPTIMIZATION OF GOLF SHAFTS
By

Nicholas Frank Abbruzzese

There are a number of tests and methods used in the golf industry that attempt to
measure and characterize golf shaft performance. The first portion of this work shows
that these tests can be modeled and studied analytically using beam finite elements. In
the second portion, the golf shaft is dynamically modeled and optimized within a finite
element simulation of a golf swing. During the optimization, parameters of the golf shaft
are varied and the performance is measured by the launch and spin characteristics of the

golf ball after impact.

To my parents, John and Terry

.0-

ACKNOWLEDGEMENTS

I’d like to thank Dr. Tom Mase and Dr. Ron Averill for their support and
guidance throughout this endeavor. I’d like to thank Dr. Mase for his patience and
understanding. His willingness and devotion to seeing me through to the end, has never
wavered. I’d like to thank Dr. Averill for his patience and unselfishness. He agreed to
take this on at much later stage in my graduate program for which I am grateful. He also
provided me with the motivation I needed at times of panic and doubt.

Most importantly, working with these two individuals has allowed me to grow in

so many ways both personally and professionally.

TABLE OF CONTENTS

 

 

 

 

 

LIST OF TABLES VII
LIST OF FIGURES I - VIII
1 INTRODUCTION AND LITERATURE REVIEW 1
1.1 INTRODUCTION ....................................................................................................... 1
1.2 LITERATURE REVIEW .............................................................................................. 2

2 SHAFT MODELING 6
2.1 OVERVIEW .............................................................................................................. 6
2.1.1 Beam Theories: Classical Euler-Bernoulli and Timoshenko .......................... 6
2.1.2 Shafi‘ EI Measurement ..................................................................................... 8
2.1.3 Assumptions of Shaft Mass ............................................................................ 10

2.2 CALCULATION OF GOLF SHAFT ZONAL FREQUENCIES .......................................... 12
2.2.1 Background ................................................................................................... 12
2.2.2 The Finite Element Model ............................................................................. 13
2.2.3 Numerical Results and Discussion ................................................................ 1 7

2.3 MODELING OF KICKPOINT .................................................................................... 21
2.3.1 Background ................................................................................................... 21
2.3.2 Buckling ~ Finite Element Model ................................................................. 22
2.3.3 Buckling ~ Numerical Results and Discussion ............................................. 25
2.3.4 Static Test ~ Finite Element Model ............................................................... 26
2.3.5 Static Test ~ Numerical Results and Discussion ........................................... 27

2.4 SUMMARY OF SHAFT MODELING .......................................................................... 29

3 OPTIMIZATION USING COMIVIERCIAL CODES 30
3.1 OVERVIEW ............................................................................................................ 30
3.2 THE “SILICON BYRON” MODEL ............................................................................ 32
3.2.1 The Shaft ....................................................................................................... 32
3.2.2 The Club head ............................................................................................... 34
3.2.3 The “Arms”, “Wrists”, and “Hands” .......................................................... 34

 

3.2.4 The Swing: Loading, Joint Definitions, Contact Definitions ........................ 35
3.2.5 The ball ......................................................................................................... 37
3.3 I-IEEDS OPTIMIZATION OF SHAFT STIFFNESS ....................................................... 38
3.3.1 Background/Aim ........................................... 38
3.3.2 Baseline Design ............................................................................................ 38
3.3.3 Preparations for Optimization ...................................................................... 43
3.3.4 Optimization Problem Definition .................................................................. 45
3.3.5 Model Variables, Responses, and Constraints ............................................. 47
3.3.6 Results and Evaluation of Optimal Designs ................................................. 55

4 CONCLUSIONS AND DISCUSSION 59
4.1 SHAFT MODELING ................................................................................................ 59
4.2 OPTIMIZATION ...................................................................................................... 59
4.3 FUTURE WORK ..................................................................................................... 60

BIBLIOGRAPHY . 61

vi

LIST OF TABLES

TABLE I. GRAPHITE DRIVER BALANCE POINTS FOR DIFFERING MASS SHAFTS. .............................................. 10
TABLE 2. MASS PER UNIT LENGTH OF THE 12 SHAFTS MODELED. ............. ' .................................................... I I
TABLE 3. STANDARD NUMBER OF E] MEASUREMENTS MADE ALONG THE LENGTH OF A GOLF SHAFT. .......... 16

TABLE 4. CRITICAL BUCKLING LOADS, LOCATIONS OF MAXIMUM DEFLECTION, AND LABELED STIFFNESS FOR

EACH SHAFT (MEASURED FROM TIP). ..................................................................................................... 26
TABLE 5. RESULTS OF STATIC TEST (MEASURED FROM TIP). ........................................................................ 28
TABLE 6. ELEMENTS AND CORRESPONDING SECTIONS OF THE SHAFT. .......................................................... 33
TABLE 7. CLUB HEAD COMPONENTS. ............................................. - - - ........................ 34

 

TABLE 8. COORDINATES NEEDED TO DEFINE POLYNOMIALS WITH RANGES USED IN OPTIMIZATION. THE X
COORDINATES RANGE OVER THE LENGTH OF SHAFT BEING OPTIMIZED (0-282 IN) WHILE THE Y

COORDINATES RANGE OVER THE RADIUS AND THICKNESS VALUES RESPECTIVELY (ALL UNITS IN INCHES).

.............................................................................................................................................................. 51
TABLE 9. CONSTRAINTS AND OBJECT IVE OF OPTIMIZATION .......................................................................... 54
TABLE 10. BALL FLIGHT AND SPIN CHARACTERISTICS OF OPTIMAL DESIGNS. ................................ 58

vii

LIST OF FIGURES

FIGURE 1. FIRST TRANSVERSE BENDING NATURAL FREQUENCY OF A SIMPLY-SUPPORTED BEAM AS A

FUNCTION OF BEAM-LENGTH ASPECT RATIO (H/L) (FROM FRIEDMAN AND KOSMA'I‘KA [7]) .................... 7
FIGURE 2. DIAGRAM OF FLEX TEST APPLIED To A SEGMENT OF GOLF SHAFT. ................................................. 9
FIGURE 3. DEVICE FOR MEASURING El PROFILE OF GOLF SHAFTS. ................................................................ 10
FIGURE 4. COMMON FREQUENCY MACHINE (FUJIKURA) USED FOR MEASURING ZONE FREQUENCIES. .......... 12
FIGURE 5. BEAM ELEMENT DEGREES OF FREEDOM. ...................................................................................... 13
FIGURE 6. SUMMARY OF MODEL DETAILS AND BOUNDARY CONDITIONS. ..................................................... 18

FIGURE 7. RESULTS FROM THE ORIGINAL CLAMPED MODEL COMPARED WITH EXPERIMENTALLY MEASURED

VALUES FOR THE 12 SHAFT S LISTED IN TABLE 2. THE ZONES CORRESPOND TO THE DIFFERENT

CANTILEVERED LENGTHS WHILE THE BARS ARE THE FREQUENCY VALUES. ........................................... 18
FIGURE 8. MEASURED CLAMP STIFFNESS OF FUJIKURA FREQUENCY MACHINE. ............................................ 19
FIGURE 9. ZONAL FREQUENCIES CORRECTED BY ADDmON OF ROTATIONAL SPRING. ................................... 20
FIGURE 10. COMPARISON OFFE MODEL SOLUTION WITH THAT FOUND BY EULER. ............................... 25
FIGURE 1 l. INmAL STARTING POINT OF SWING SIMULATION. ...................................................................... 39
FIGURE 12. SIMULATION AT 0.15 SECONDS. ................................................................................................. 39
FIGURE 13. SIMULATION AT 0.30 SECONDS. ................................................................................................. 40
FIGURE 14. SIMULATION AT 0.37 SECONDS. ................................................................................................. 40
FIGURE 15. SIMULATION AT 0.41 SECONDS. ................................................................................................. 41
FIGURE 16. SIMULATION AT 0.42 SECONDS. ................................................................................................. 41
FIGURE 17 . SIMULATION JUST AFTER IMPACT AT 0.43 SECONDS ................................................................... 42
FIGURE 18. SIMULATION AT 0.44 SECONDS. ................................................................................................. 42
FIGURE 19. STIFFNESS PROFILE OF BASELINE SHAFT DESIGN. ....................................................................... 43
FIGURE 20. SEQUENCE OF SWING SIMULATION. ............................................................................................ 45
FIGURE 21. VARIABLES INFLUENCING BACKSPIN .......................................................................................... 46

viii

FIGURE 22. AN EXAMPLE POLYNOMIAL FOR THE OUTER RADIUS OF THE SHAFT. THE MAXIMUM AND
MINIMUM ALLOWABLE RADIUS IS ALSO SHOWN. ................................................................................... 49
FIGURE 23. AN EXAMPLE POLYNOMIAL FOR THE THICKNESS OF THE SHAFT. THE MAXIMUM AND MINIMUM
ALLOWABLE THICKNESS IS ALSO SHOWN. .......................................... ' .................................................... 50
FIGURE 24. SEQUENCE OF ONE EVALUATION IN OPTIMIZATION. ................................................................... 53
FIGURE 25. OPTIMAL DESIGNS - GROUP 1. THE PROFILE IS THAT OF THE 28.2 INCH SPAN VARIED IN THE

OPTIMIZATION. ...................................... .. - - ......................... - - - - .................. 56

 

 

FIGURE 26. OPTIMAL DESIGNS - GROUP 2. THE PROFILE IS THAT OF THE 28.2 INCH SPAN VARIED IN

OPTIMIZATION. ................................................................................................................... 57

 

FIGURE 27. DESIGN 65 EXHIBITING A HIGH DEGREE OF DEFORMATION IN DOWNSWING. .............................. 58

1 Introduction and Literature Review

1.1 Introduction

In recent years, the golf industry has been searching for more objective and
technical means of characterizing golf shaft performance. Conventional terms and
methods used in the industry have been shown to not accurately describe shaft
performance; they remain too subjective and the methods by which they’re determined
lack consistency throughout the industry [1]. AS a result, attempts are being made to
better quantify and “fit” the performance of a certain shaft to a specific golfers swing.

Currently there are two basic attributes used in evaluating golf shaft performance.
These two attributes are also the cornerstones on which the shafts are marketed and sold.
They consist of measures of the shafts flex ,or“stiffness,” and “kickpoint”. There are a
wide range of experimental tests and methods used to measure these characteristics.
What Chapter 2 will Show is that a majority of these experimental tests can be accurately
modeled using elementary beam finite element models. This is accomplished through an
initial knowledge of the shafts “bending profile” or “E1 curve”. The “bending profile” or
“E1 curve” is becoming a standard in the industry for benchmarking shaft designs [2].
Prior knowledge of this bending profile is what the elementary beam models rely upon.
This will be discussed in further detail in Chapter 2.

Stemming from the bending profile and its increasing importance to the overall
evaluation and design of golf shafts, an optimization study is also undertaken to find the
optimal bending profile for a given golf swing. Certain design variables of the golf shaft

are dynamically studied and optimized using a full—swing finite element simulation in LS-

DYNA. Unlike the static analyses of Chapter 2, the complete dynamic response Of the
shaft is evaluated and optimized. The details of the finite element model and of the
optimization are presented in Chapter 3. The following section will provide a broad

technical review of golf shafts.
1.2 Literature Review

Perhaps the best paper to lead this review is by GP Horwood, “Golf Shafts — a
technical perspective” [3]. Horwood acknowledges the impact of the golf shaft on the
overall golf swings performance, while also acknowledging the extreme difficulty with
finding consistent and technically sound ways of characterizing this performance. He
warns of the sometimes extreme marketing hype that is used in the industry which leads
golfers into believing that there are certain shafts fitted “perfectly” for their individual
swings. He explains on a very simple level, that the characteristics of a golf shaft that are
of most importance are the golf shaft’s bending and torsional stiffness- In the golf swing
the inertia of the club head coupled with the offset of the club head center of gravity are
what influence the golf shaft dynamics in the golf swing.

In the paper by Chou and Roberts the accuracy of three popular methods for the
measurement of a golf shafts “kickpoint” are analyzed [1]. In their Study three steel and
three graphite golf shafts with a Wide range of kickpoints were chosen. The six drivers
were then subjected to machine testing. The goal was to relate the changes in ball
trajectory to the variations in shaft kickpoint. Classical golf theory has always suggested
that the degree and location of a golf shafts kickpoint has a direct influence on ball
trajectory. However, the conclusions made by Chou and Roberts showed a poor

correlation.

Wallace and Hubbell examined the influence of golf shaft stiffness on the swing
kinematics and dynamics of a group of amateur golfers [4]. Each golfer was asked to hit
3-10 Shots with each of 3 different S-irons varying only in shaft stiffness. Variables
considered in the study included club head speed, handicap, ball speed, clubface angle,
swing-path angle, as well as a few others. Although this study provided and reaffirmed
some of the general conclusions made about the differences between low and high
handicappers, it failed in drawing any sound conclusions in regards to shaft stiffness.
This was due to analyzing the golfers as a group, rather than looking at how Shaft
stiffness might affect an individual golfer. Wallace and Hubbell also acknowledged that
the study did not take into account important launch conditions like ball launch angle and
spin rate.

The work by Butler and Winfield attempts to find a better understanding of shaft
dynamics in the downswing along with finding key variables that affect the dynamic
performance [5]. Their study consisted of a finite element model along with certain
experimental measurements of shaft deflection and twisting during the downswing. The
finite element beam model developed by Butler and Winfield helped to Show the
influence of the club head’s Offset center of gravity to the twisting of the shaft. The first
mode of the shaft was analyzed and a simple correlation was made between the deflection
of the shaft and the coupled twisting of the club head. The most significant portion of the
paper was the experimental measurements that provided some very basic insight into the
shaft deflection in the downswing. Measurements of the “lead/lag” deflection which
occurs in the actual plane of the golf swing and measurements of the “toe up/down”

deflection occurring in the out of plane direction were made for 3 varying swing types.

The conclusions certainly helped in finding a better understanding of the dynamics of the
golf shaft. Butler and Winfield also acknowledged the importance and possible use of
such measurements for club-fitting.

Brouillette details three non-destructive methods of measuring the flexural
rigidity distribution along the length of a golf shaft [6]. These are alternative and out-
dated methods for measuring the El profile referred to in Chapter 2. The methods
discussed by Brouillette are sound and may have influenced the current method used in
the industry today.

Mase documents a new method, which is quickly becoming industry standard, for
measuring the E1 profile of golf shafts [2]. The paper details a method for correction of
the radial deflections when measuring the bending stiffness profile or E1 profile of golf
shafts (see Sec. 2.1.2).

The work by Friedman and Kosmatka is mentioned here to support the elementary
beam modeling in Chapter 2 [7]. The models were developed using simple Euler-
Bemoulli beam finite elements. Initially the frequencies of the shorter beam models,
those having a somewhat larger aspect ratio, were found to be too stiff. It was originally
thought that due to the beams shorter length that the shear effects could not be neglected.
This led to the investigation of Timoshenko beam finite elements. Figure l in Chapter 2
is referenced from the study. This figure proves that even the shortest shaft lengths
examined in the study can be accurately modeled by Euler-Bernoulli finite elements.

The work by Lee and Kim is one of the very few Optimization studies directly
applied to golf equipment [8]. In it they detail an optimization method for graphite golf

shafts. The objective function is to minimize the torque of the shaft, while satisfying

certain constraints like flex, “kick point”, and weight requirements. This is accomplished
using a finite element model developed in ABAQUS interfaced with genetic algorithm

type optimization software specifically for composite materials called DARWIN.

2 Shaft Modeling

2.1 Overview

There are a number of technical tests and methods used in the golf industry that
attempt to measure golf shaft performance. Bending stiffness, torsional stiffness,
kickpoint, durability, and conformity to the Rules of Golf are major factors in shaft
design and manufacture. In this thesis, two main areas of shaft performance were
considered: a measure of the shaft’s overall stiffness behavior, and a measure of the
shafts “kickpoint” or “flexpoint” behavior.

What this chapter will Show is that the majority of these tests, which in some
cases involve rather intensive and elaborate experimental setups, can be accurately
modeled analytically by using beam finite elements. Conclusions drawn from these
models can also be used to validate or reject their relevance in describing shaft behavior.
And in Some cases the models can even eliminate the need for experimental testing,
thereby reducing the time and cost of such testing.

The following sections provide an introduction to the models. Basic beam theory
is reviewed and important assumptions and key points in the modeling are discussed.
The models were implemented using the MATLAB program version 7.0.0.19920 (R14).
2.1.1 Beam Theories: Classical Euler-Bernoulli and Timoshenko

The Euler-Bernoulli and Timoshenko beam theories are two of the most
commonly used beam theories in structural mechanics problems. However the

underlying assumptions found in the two theories are slightly different.

In the Euler-Bernoulli theory shear effects are neglected and as a direct result the
cross-section of the deformed beam must remain plane and normal to the longitudinal
axis. On the other hand, the Timoshenko theory includes shear effects, therefore
allowing the cross-section to remain plane but not necessarily normal to the deformed
axis. For slender beams, those defined as having a beam-thickness-to-length aspect ratio
much less than one-tenth (h/L<O.l), the Euler-Bernoulli beam theory tends to give good
approximations because transverse shear deformation is negligible. However, as the
beam becomes thick and non-Slender (O.l<h/L<l.0) transverse shear deformation must be
included and the Timoshenko model should be considered [7]. The figure below is

borrowed from Reference [7] and summarizes their conclusions.

 

 

 

 

 

 

 

 

10 : ' V ' Y f T r*r—r v v v r 7 Y Y 1 I a f
-- — Bemoulli-Euler //’
- Timoshenko Soln / ]
I /
0 / W ”LL. __
1 L _,-----~— I
0 I / I“, .--"’ — :
I /,/// -1
/ 4
3 /
3
L
.1
10 T 1
v"
10.2 I I I I I I I 41 I I L 1 1 I 4 1 1
.2 -1 O
10 10 10
h/L

Figure 1. First transverse bending natural frequency of a simply-supported beam
as a function of beam-length aspect ratio (h/L) (From Friedman and Kosmatka [7])

At the outset of the study there was doubt as to which theory was most
appropriate for modeling the shorter golf shaft lengths. The initial results from the
models of the Shorter shaft lengths were overly stiff and it was thought that the
Timoshenko theory would provide better results. However, according to the figure
above, it was concluded that the traditional Euler-Bernoulli theory should be sufficient
and accurate. The overly stiff predictions were later attributed to the boundary
conditions, as discussed in the sequel.

When considering the cantilevered golf shafts in this study, the shortest length is
15 inches while the shaft diameter is approximately 0.5 inches. This gives a beam-
thickness-to-length aspect ratio no larger than 0.03. Therefore, the Euler-Bernoulli beam
theory is assumed appropriate and is used in the development of all finite element
models.

For a more detailed look at the finite element formulations of both the classical
Euler-Bernoulli theory and the Timoshenko theory, see Reddy [9].

2.1.2 Shaft EI Measurement

A fairly technical approach to help characterize shaft performance is to measure
the shaft’s “bending profile” or “E1 curve” (E = modulus of elasticity, I = moment of
inertia or second moment of area). As will be seen, this EI curve is the key to the
development of all finite element models in this part of the study. The data found from
these measurements are direct inputs to the finite element models. A brief description of
the method used in obtaining these measurements is now discussed.

The modulus of elasticity is a measure of material stiffness and moment of inertia

is a measure of material distributiOn. Together, they make up a structural Stiffness per

unit length. True stiffness, EI/L, is not used Since a constant span length of 250 mm is
used in the experimental apparatus. The profile consists of a series of bending
measurements taken along the length of the shaft. From end to end, incremental spans of
the shaft are subjected to what is called a three point flex test; two ends are simply
supported while a point load is applied at the mid-span. Under this type of loading, a
Simple mechanics of materials relation can be used to provide a local estimation of the

shaft’s bending stiffness or E1.

_ PL2

EI -
486,

 

(2.1-1)

Here P is the applied load, L is the Span length, and 5c the beam deflection at the mid-
span. Figure 2 is an example showing a schematic of the flex test. The actual test

machine for measuring shafts is Shown in Figure 3.

SON SON

 

‘2‘

“A

.1... -L A m . u win:
I‘— 125 m J: 125 mm ——>I

lOON

 

Figure 2. Diagram of flex test applied to a segment of golf shaft.

 

Figure 3. Device for measuring EI profile of golf shafts.

These bending profiles provide insight into how the shaft’s bending stiffness changes

along its length. This information becomes the foundation for modeling the golf shaft.

2.1.3 Assumptions of Shaft Mass

Another key assumption made in the modeling is in regards to the golf shaft mass.

The mass is assumed to be distributed evenly along the length of the shaft. Table 1

shows the relationship between the shaft mass and its distribution about its balancing

point. These shafts were selected in order to include the range of weights and designs

currently available. The numbers reported are an average of n different shafts that were

measured.

Table 1. Graphite driver balance points for differing mass shafts.

 

Shaft

Mass, 1 Balance Point %

 

Graphite Shaft 1 (n=4)
Graphite Shaft 2 (n=4)
Graphite Shaft (n=3)
Graphite Shaft (n=3)

$175 51.9
65.5 52.2
70.3 51.4

32 53.1

 

The mass of both steel and graphite Shafts are distributed quite evenly about the
balance point. Further support for this assumption stems from the inverse relationship
between a typical shaft diameter and its thickness.

In the subsequent finite element models, the mass per unit length pL is calculated
as an average for each shaft and is found by dividing the recorded shaft mass by the
overall length. Table 2 shows the 12 shafts used in the models and their assumed mass

per unit length.

Table 2. Mass per unit length of the 12 shafts modeled.

 

Length Weight Mass per unit

 

M°°°' F'” (In) (9) Length (kg/m)
YS-6 Fl 46 63 0.0539
YS-6 S 46 64 0.0548
YS-6 SX 46 65 0.0556
YS-6 X 46 65 0.0556
YS-6 TX 46 68 0.0582
YS-7 H 46 76 0.0650 '
YS-7 S 46 77 0.0659
YS-7 3X 46 78 0.0668
YS-7 X 46 78 0.0668
YS-7 TX 46 79 0.0676

YS-8. 1 S 46 86 0.0736

YS-8.1 SX 46 87 0.0745

11

2.2 Calculation of Golf Shaft Zonal Frequencies

2.2.1 Background

The shaft’s first natural frequency at various cantilevered lengths are measured
and are referred to as “zonal frequencies” in industry. The frequencies at cantilevered
lengths of 991, 787, 584, and 381 mm (39, 31, 23, and 15 inches respectively) are
determined experimentally using a modified Fujikura frequency machine (Figure 4).
These “zonal frequencies” are being used as yet another way to characterize shaft

stiffness and performance.

 
    
 
 

Light Sensors Clamping
Air Cylinders

  

205 g Weight

Figure 4. Common frequency machine (Fujikura) used for measuring zone
frequencies.

Currently, as the popularity of these new methods increases and as new shaft
designs are constantly emerging, large databases must be kept to include the measured
data. The following willshow that a simple beam finite element model can accurately

predict the zonal frequencies from the measured bending profile. Therefore reducing the

time and cost associated with experimentally measuring the zonal frequencies. The finite
element formulation and complete description of the model will be highlighted followed
by a comparison of numerical results for 12 shafts.

2.2.2 The Finite Element Model

The golf shaft is modeled as a cantilever beam with a point mass included at the
tip to represent the club head. Like most cantilever beams, the original model was
clamped at one end having both the transverse and rotational displacement of one node
constrained. However, as it will be explained in the results, this original model had to be
modified. The constrained end of the shaft was remodeled to include a simple support
along with a stiff rotational spring. The simple support continued to constrain the
transverse displacement, but the rotational spring was needed to model a small amount of
compliance in the rotation of the joint.

The discretization of the beam is accomplished using a number of Euler-Bernoulli
finite elements. Use of simple Euler-Bernoulli beam elements was justified in Section
2.1.1. A typical depiction of the element can be seen in Figure 5. Each element has two
nodes each with two degrees of freedom; a transverse displacement w, and a cross-section

rotation 0.

6' W1 W2 02

Cl; . ___:|+l>

Figure 5. Beam element degrees of freedom.

 

 

 

l3

The formulation of the stiffness matrix and the mass matrix for a single element
can be found in many introductory finite element textbooks [9,10]. It will be briefly
summarized here for the sake of completeness.

The strong form of the Euler-Bernoulli free vibration equation is given as

2 2W
p 3”’+ a —22[Ela w]: o (2.2-1)

 

 

L 812 +x28 Ox
where pl, is the mass per unit length, E1 the bending stiffness, and w(x,t) the transverse

displacement of the beam. The strong form here requires the solution of w(x,t) to be of

the C4 type having continuous 4th order derivatives. The weak form, as it will be shown,
will allow this restriction to be “weakened” to require a solution of C2 type, which
requires only continuous 2“d order derivatives.

To arrive at the weak form of the equation, Equation 2.2-l is multiplied by an

arbitrary test function Fix, t) and integrated over the domain. Provided that the test

function is zero wherever an essential boundary condition exists; the final result

simplifies to

32w 02w
28x2

 

 

jag; tywdx+IEIax2 dx=0 (2.2-2)

where the domain is considered over the entire length of the member.
The next step is to then discretize the weak form of the equation. The

displacement across each element is assumed to take on the form

w(x,t) = N, (x)w, (t) + N,(x)r9I (t) + N3(x)w2 (t) + N, (no, (t) (2.2.3)

14

where N1, N2, N3, and M are Hermitian cubic interpolation functions (also called shape
functions). These are time independent and are used to ensure continuity across the
elements.

Using the discretization above in the weak form yields the finite element mass

and stiffness matrices. The entries are given by the following

 

L
M q, = IpLNa(s)Nfl(s)ds a,,B =1,2,3,4 (2.24)
0
L 2 OZN
KMap = IEIa 322(3) 352(3) 618 61,,6 = I,2,3,4 (2.2-5)

0

H 99
S

where a local coordinate system has been introduced for ease in computing the
relative entries.

The resulting consistent mass and material stiffness matrices for a single Euler-

Bernoulli finite element are

 

 

' 12 6L —12 6L '
E1 4].2 —6L 2L2
K =— 2.2-6
M L3 12 -6L ( )
Lsymmetric 4L2 d
and
’ 156 22L 54 —I3L'
4L2 13L -3L2
= PLL (2.27)
420 156 —22L
Lsymmetric 4L2 3

 

 

where L is the element length, and the mass per unit length p1. and bending stiffness E]
are assumed constant over the element.
AS seen in the matrices above, approximations in the model are needed for each

element’s bending stiffness E1, and mass per unit length p1,. The latter of the two is

15

 

approximated by dividing the total mass of the Shaft by its length, giving an average
value of p], for all elements. This is a fair approximation because the mass along the
length of the shaft can be considered relatively constant (seeSec. 2.1.3).

The approximation for each element’s bending stiffness comes from prior
knowledge of the shafts bending profile (see Sec. 2.1.2, [2]). This measured bending
profile or E1 curve provides approximate values of the shaft’s bending stiffness along its
length. The shaft is flexed over a 250 mm (10 in) span and a load ranging from 67 to 100
Newtons (15-23 lbs) is applied to the mid-span. This procedure is repeated every 50 mm
(2 in) along the shaft up to 975 mm (39 in). The location of the point load is taken as the
point at which the bending stiffness is evaluated. The table below details the location of

each measurement taken along the shaft.

Table 3. Standard number of E1 measurements made along the length of a golf

 

 

shaft.
Location of
El Measurement Applied Region
Measurement from Tip End of Shaft in
Number (mm) Models (mm)
1 125 0 - 150
2 175 150 - 200
3 225 200 - 250
4 275 250 - 300
5 325 300 - 350
6 375 350 - 400
7 425 400 - 450
8 475 450 - 500
9 525 500 - 550
10 575 550 - 600
11 625 600 - 650
12 675 650 - 700
13 725 700 - 750
14 775 750 - 800
15 825 800 - 850
16 875 850 - 900
17 925 900 - 950
1 8 975 950 - 1000

 

16

Also shown in Table 3 are the regions of shaft in the model over which each
measurement is assumed to apply and be constant. For example, the first measurement
taken near the tip of the shaft at 125 mm is applied to the first 150 mm of the shaft in the
model. These regions helped determine where elements should be divided and placed.
AS it turned out, dividing each region into two elements was sufficient enough for the
frequencies to converge. All elements, except for those at the very ends of the model, are
25 mm in length.

Various cases were studied to test the validity of the models meshing as well as
the logic used in assigning the E1 measurements to regions. The original mesh was
refined to contain twice the number of elements and there was a negligible change in the
frequencies. In another case, the El measurements were fit with a 4th-degree polynomial
aiming to smooth the transition of the E1 values between elements. This was also shown
to have a negligible change in the resulting frequencies.

2.2.3 Numerical Results and Discussion

Given the bending profile data for 12 graphite-epoxy Shafts, the zonal frequencies
for each Shaft could be computed using the finite element model presented. Results from
the original clamped model and results from the model using the rotational spring are

presented. Details of the various models are shown in Figure 6.

l7

Different Boundary

Conditions MOdCI Details

 

label-991mm

36 beam elements
b M
-——§ I

Cantilevcred Zone 2 - 787 mm

end 28 beam elements i
Zone 3 — 584 mm
=@ 20beemelements i

Simple end
with torsion Zone 4 - 381mm

spring 18 beam elements Q

Figure 6. Summary of model details and boundary conditions.

 

 

The results from the original clamped model can be seen in Figure 7. The

computed zonal frequencies are compared to the experimentally measured frequencies.

 

 

Computed Zonal Frequencies Measured Zonal Frequencies
. . . r 9]] . . . .
133515
800 - 81]) - .
700 - 7CD ~
600 - 5m .
’E‘ ’E‘
D. D.
8 500 ~ 3 am -
3 F
5 5
as; 400 g 40.1 -
u‘: u‘:

EB
88

.-
D
O
_n
8

 

 

 

 

   

Figure 7. Results from the original clamped model compared with experimentally
measured values for the 12 shafts listed in Table 2. The zones correspond to the
different cantilevered lengths while the bars are the frequency values.

18

The computed frequencies become overly stiff as the shaft length shortens. The
average error for the first zone is less than 1 percent, while the error for the second, third,
and fourth zones increases to 2.44, 5.94, and 13.35 percent, respectively. This led to a
closer inspection of the experimental equipment used to measure the frequencies.

The Fujikura machine was found to have a small amount of compliance in the
clamped joint. The joint was in-effect, only semi-rigid. The air cylinders used for
clamping the shaft make measurements easily and quickly attainable, but build in this
small amount of compliance. Therefore, the original finite element model was modified
to account for this small, rotational flexibility and a rotational spring was added to the
end node. The stiffness of the clamp was measured and the following figure shows the

linear fit of the stiffness value.

 

 

 

 

Fujikura Clamp Stiffness
12 I I I I fl I f j T
/6
. / _
10 /
/
d
8 / I
//
E
' 9’
E 5 x/
I: - -I
(D
/
E x
E ’0'
4 ~ / -
//
o/ Itm = 127349-1004 N-m
/
2 P / -
/
/
/
l l l l I l l l l
08 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Rotation, rad x 10-3

Figure 8. Measured clamp stiffness of Fujikura frequency machine.

I9

Using the measured rotational stiffness in the finite element models greatly
improved the error in the measurements of the shorter lengths without altering the
accuracy of the longer lengths. Figure 9 Shows the corrected frequency plots with the

average error listed above each zone.

 

 

 

 

 

 

 

 

Computed Zonal Frequencies Measured Zonal Frequencies
BOO - 7.35% . 3m . ‘
700- H — 700- -
500- i - see» I] _
’5 I ’E‘ 1
9. ll 9. .
3500 3.87% In . 3500- 3‘ _
>“ [I >. I
0 1 ‘I, U
5 II I] g _
3 400 ~ ' I] . 3 4m . II}. _
S 1.35% i g" l‘r ,
I: , , I I: ‘l‘
300 0.39% [‘H II ‘: - 3m . i H ‘ _
"I l I ll II I” ‘ ‘ 1
‘l " I! . ‘ .
200‘ I 11; I =l ‘ 200‘ ”I i l» w -
I . l} I II I ‘ l} 31 '
i I I l .‘l ‘i . I! ‘
100» ‘i. [I I ll - 1013- ll‘ , I] . 4
5‘; ‘ ‘l ‘ I ‘ I ::3 ‘
Ti; Ii l H ' i l I], ‘
1

 

2
Zones Zones

Figure 9. Zonal frequencies corrected by addition of rotational spring.

Although the error at the shorter shaft lengths does improve significantly, it still
falls short of accurately predicting the frequencies for zones 3 and 4. Further modeling
and investigation of the clamping mechanism is needed. Other sources of error
contributing to this discrepancy could involve the experimental measurement device.
The light sensor may have difficulty measuring the higher frequencies, therefore under-

prcdicting the measured results.

20

2.3 Modeling of Kickpoint

2.3.1 Background

A golf shaft’s “kickpoint”, also referred to as “flex point” or “bend point”,
remains an issue of much debate in the golf industry. There are many terms, and just as
many methods that are used in hopes of characterizing this aspect of golf shaft
performance. The goal is to effectively relate the degree and location of a golf shaft’s
“flex” during the downswing to the resulting ball flight and trajectory. It is thought that a
shaft with high “kickpoint” will result in lower ball flight, whereas a low “kickpoint” will
result in higher ball flight.

One of the many proposed methods for determining a shaft’s “kickpoint” is by
analyzing its buckling behavior. This test consists of applying an increasing axial load to
the shaft until buckling is reached. The evaluation of the shaft’s “kickpoint” is then
found from the buckled shape of the shaft. The maximum point of deflection and the
corresponding location are then used to quantify and characterize the shaft “kickpoint”.

Yet another method that will be modeled in the sections that follow consists of a
simple static test. The shaft is clamped near the butt end and a specified load (usually of
2-lb or 8.9N) is hung from the tip of the shaft. A string is then connected between the
two ends of the shaft creating a line of reference. The point at which the shaft is furthest
from the string is noted and used as a relative evaluation of the “kickpoint”. The static
test remains controversial and has fallen short of being an accurate means of describing

“kickpoint” as the finite element model will show.

21

2.3.2 Buckling ~ Finite Element Model

By definition, buckling is the bifurcation or branching of equilibrium states from
stable to unstable. The analysis most commonly used in engineering practice with
respect to buckling is what is known as linearized prebuckling analysis (LPB) [11,12]. In
this analysis, solving a traditional linear algebraic eigenvalue problem involving the
tangent stiffness matrix K. can assess the stability of a structure. The resulting
eigenproblem yields the magnitude of the load at which buckling occurs (lowest
eigenvalue) as well as the corresponding buckled shape of the column (corresponding
eigenvector). The key steps and assumptions in LPB will be briefly summarized.

In LPB analyzing the tangent stiffness matrix K. assesses the stability of the
structure. Here p.- and z.- are the im eigenvalue and eigenvector, respectively, of the
matrix.

K,z,. = I,u,.z,. ‘ (2.3-l)
Since the tangent stiffness matrix is real symmetric all of its eigenvalues are therefore
real. However as certain parameters are varied, in most cases the loading of the
structure, the eigenvalues ,u,- of the tangent stiffness matrix change and the matrix. can
become singular.

K‘ = K,(}l) (2.3-2)
It is at this point, where the system transitions from a region of strong stability (p.- > 0) to
neutral stability and singularity (p,- 2 O) and into instability (u,- < O). The point at which

the matrix becomes singular is therefore considered to be the critical value denoted as A".
det K, (1”) = O (2.3-3)

Loading the structure beyond this value initiates bifurcation and the structure will buckle.

22

The exact formulation of the buckling problem involves the decomposition of the
tangent stiffness matrix K into the material and geometric stiffness matrices KM and Kc,
respectively.

K =KM +KG (2.3-4)
The material stiffness KM is a function of the physical properties in the element and is
constant, whereas the geometric stiffness K0 is dependent upon the load in the element.
As noted above the stability test requires that K, be singular, which leads to the final form

of the stability eigenproblem.
K,z,. =(KM +ZiKG)zi =0 (2.3-5)

The critical loads required for buckling can now be obtained directly from the solution of
the above equation. The eigenvalue 1; closest to zero is the lowest critical load needed to
reach buckling, while the corresponding eigenvector z,- is the buckling mode shape.

In the analytical model presented, each shaft is again discretized into a sufficient
number of classical Euler-Bernoulli beam elements. The shaft length used in the model
is 1.05m (41.5 in). The element material stiffness matrix was given previously in
Equation 2.2-6 but is repeated here along with the exactly integrated geometric stiffness

matrix [10,12].

23

 

 

 

" 12 6L —12 6L7
4L2 —6L 2L2
KM = £31 12 _6L (2.3-6)
Lsymmetric 4L2_
and
' 36 3L —36 3L“
4L2 —3L —L2
KG Li. (2.3-7)
30L 36 —3L
Lsymmetric 4L2_

 

where P is the axial force in the element. In the simple case of modeling the golf shaft,

the structure is statically determinate therefore allowing P to equal the applied load.

With the element matrices defined, the final steps involve assembling the global
matrices and applying constraints to the model. The ends of the shaft are constrained by

simple supports or pins eliminating their translational degrees of freedom.

The accuracy of the model was verified through comparison with the well known

Euler’s method, published in 1744. For a uniform beam/column of constant material and

geometric properties subjected to an axial load, Euler’s Formula gives the lowest critical

buckling load.
(2.3-8)

 

The beam used as a test case in the model had a bending stiffness of E1 = 50 Nm and a

length of 1.05m (41.5 in).

24

buckled mode shape

 

 

 

 

 

 

 

 

0.5 7 T I I I I T T I
0.45 ~ /,»\ ~
/ \ .
/ \
o 4 — / \ .
/ \
‘5 / \
‘1’ 0.35 ~ / , ~
E - \
Q)
8 /
E- 0.3 " / "l
U /
§ 0.25 ~ / —
.C
U)
§ 0.2 ~ / _
E
0 15 / \‘ —
' / .- FE solution \
0 1 j/ — Theoretical w=A'sin(n*pi'x/L) \\
\
0.05 I l J l 1 l 1 1 l J
O 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
MGIGI'S

Figure 10. Comparison of FE model solution with that found by Euler.

Figure 10 shows the buckled mode shape compared with the theoretical solution
obtained by Euler. The comparison of the critical loads is shown below where the
superscript E represents the solution found using Euler’s Formula.

F" = 444.126N Ff = 444.126N (2.3-9)

The results from the finite element model are in complete agreement with the
results found by using Euler’s Method.
2.3.3 Buckling ~ Numerical Results and Discussion

Given the bending profile data for the group of 12 shafts, the critical buckling

loads Fe, and the location of maximum deflection could be computed using the finite

element model and LPB analysis presented.

25

Table 4 shows the calculated critical loads for each shaft, the location of
maximum deflection from the tip end, as well as each shafts industry labeled stiffness.
As seen, the magnitude of the critical load is directly related to the stiffness of the shaft

but there is no discemable relation with the points of maximum deflection.

Table 4. Critical buckling loads, locations of maximum deflection, and labeled
stiffness for each shaft (measured from tip).

 

 

Pt of Max.

Deflection Stiffness
Fcr (N) (m) Label
279.75 0.456 R
305.43 0.459 S
316.66 0.448 SX
343.1 1 0.453 X
365.18 0.456 TX
274.55 0.464 F1
31 1.89 0.464 S
330.91 0.464 SX
332.57 0.465 X
339.69 0.468 TX
302.58 0.46 S
325.28 0.462 SX

 

The range of values shown for the locations of maximum deflection is very small,
therefore making the characterization of the shafts “kickpoint” very difficult. This fact is
well known to those in the golf industry.

2.3.4 Static Test ~ Finite Element Model

The static test for determining flex point can be modeled as a simple cantilever
beam subjected to a tip load. The unknown displacements w are found from the solution
to the equation below

Kw = F (2.3-10)

where K is the global material stiffness matrix and F the equivalent nodal force vector.

26

As in the other models, the shaft was discretized into a sufficient number of
Euler-Bernoulli finite elements. The element material stiffness matrix K,,. was given
previously in Equation 2.2-6, and the global matrix K was assembled from each element

stiffness matrix Km.

The only applied force in the model is the tip load of 0.907-kg (2-lb), which was
modeled as an equivalent nodal force in the transverse direction. The equivalent nodal
force vector therefore consists of only one non-zero entry corresponding to the node at
the tip of the shaft.

The accuracy of this model could be verified through comparison to a well-known
strength of materials solution. Given a uniform cantilever beam subject to a tip load, the

exact solution for the maximum displacement is given by

PL3
WMAX = _i—I-

(2.3-11)
whereP is the applied load, L the length of the beam, E the young’s modulus, and I the
moment of inertia. The beam used to validate the model had a length of L = 1.05m, load
of P = 10.0 N, and combined stiffness E1 = 50 N'mz. The maximum displacement was
calculated to be -0.0772 m and compares exactly with the result from the strength of
materials solution.
2.3.5 Static Test ~ Numerical Results and Discussion

Given the bending profile data for the group of 12 shafts, the location of

maximum deflection could be computed using the static finite element model presented.

Table 5 summarizes the results.

27

Table 5. Results of Static Test (measured from tip).

 

 

Pt of Max.
Deflection Stiffness
(m) Label
0.537 R
0.54 S

0.531 SX
0.536 X
0.539 TX
0.543 H
0.542 S
0.543 SX
0.544 X
0.547 TX
0.538 S
0.541 SX

 

As with the buckling results, the static test results are also inconclusive and lend

little insight into characterizing or making any kind of distinction regarding “kickpoint”.

28

2.4 Summary of Shaft Modeling

It has been shown that many of the golf industry tests used to characterize golf
shaft performance can be modeled and studied. With prior knowledge of the shaft’s
bending profile or E1 curve, a number of these industry standard tests can be modeled.

' The zonal frequencies test, which serves to characterize overall shaft stiffness,
was first modeled. The calculation of the zonal frequencies matched well with the
experimentally measured values at the lower shaft zones, with a slight over-prediction at
the higher zones. Further investigation of the Fujikura frequency machine is
recommended in order to pinpoint the cause of this error.

The purpose of the next two tests modeled was to help characterize shaft
kickpoint. This term is rather loosely defined and its measurement has never been
standardized across the industry. It was shown that modeling of these two tests is
possible; however their inconclusive results proved that their relevance in describing this

property of golf shafts is rather poor.

29

 

3 Optimization using Commercial Codes

3.1 Overview

In the previous Chapter it was shown that certain tests and methods used for
measuring and quantifying golf shaft performance can be easily modeled and studied
analytically. The characteristics of the shaft’s stiffness and bending properties were
studied without concern for the dynamics of the shaft in an actual golf swing. In the
current Chapter a broader study is undertaken, in which the shaft is dynamically modeled
and optimized using the full-simulation of a golf swing. During the optimization,
parameters of the golf shaft are varied and the performance is measured by the launch
and spin characteristics of the golf ball after impact.

Dr. Tom Mase of the Composite Materials & Structures Center at MSU provided
the initial finite element model used in the study. It was originally created by Dr. Mase
to model the mechanical swing machines used in the golf industry affectionately known
as “Iron Byron”. However it soon became apparent that the model could also be utilized
in an optimization study of the golf shaft bending stiffness. The initial model was then
modified accordingly.

The swing is simulated using LS-DYNA, which is a commercially available finite
element software capable of nonlinear dynamic analyses. The club and golf swing are
modeled as well as the golf ball and its impact with the club during the swing. LS-
DYNA was chosen for this simulation because of its specific capabilities of modeling
contact-impact interfaces. The software originated from a need for the stress analysis of

structures subjected to various impact loadings. The overall model consisted of the ball,

30

the club head and shaft, and rigid connections used to represent the arms, wrists, and
hands. Details of the model will be fully explained in the following sections.

The optimization of the golf shaft was accomplished. through interfacing the LS-
DYNA swing simulation with the commercial optimization software HEEDS
(Hierarchical Evolutionary Engineering Design System). Parameters of the shaft
stiffness were changed in each swing by HEEDS and each evaluation’s performance was
measured until optimal solutions were reached. Each run of the swing‘simulation was
evaluated under certain constraints and objectives concerning the golf ball flight and spin.
Further details and the results of the optimization will be discussed in later sections as

well.

31

3.2 The “Silicon Byron” Model

3.2.1 The Shaft

The golf shaft was modeled in LS-DYNA using 3-node beam elements. The
material was declared as linear elastic. There were a total of 37 beam elements, which
made up 32 distinct sections (Table 6). Each section contained specific geometric
properties that were used as variables in the optimization study. Details on the use of
these sections and variables will be discussed in detail in the following section on
optimization. The particular initial design for the shaft was taken as a long time popular
steel stiff shaft (True Temper Dynamic Gold S300). In the next subsection, the parts of
the robot are briefly described. Following that, a subsection is presented that describes
the time varying moments that were applied to the rigid bodies to drive the model and

thus swing the club.

32

Table 6. Elements and corresponding sections of the shaft.

 

w u.._—---. ..-_- 7 , .. 7.,-..2__._.1-2--_._._.._

 

Element # Section ID # of Section

1

2 4 6
3

4 3
5

5 3
7

8 * 1
9 * 1
10 * 9 1
11 " 10 1
12 " 11 1
13 * 12 1
14 * 13 1
15 ’ 14 1
16 ' 15 1
17 * 16 1
18 ’ 17 1
19 ' 18 1
20 ' 19 1
21 * 20 1
22 * 21 1
23 ’ 22 1
24 " 23 1
25 * 24 1
26 * 25 1
27 * 26 1
28 * 27 1
29 ’ 28 1
30 * 29 1
31 * 30 1
32 * 31 1
33 ' 32 1
34 ’ 33 1
35 * 34 1.2
36

37 35 3.8

* Elements used as variables in the. antimigatiQUMM

 

 

33

3.2.2 The Club head

The club head was modeled after the majority of drivers currently in the industry
(circa 2004) having a volume of 360 cc. It consisted of eight different parts all of which
were meshed using 4-node shell elements. Different parts were used to model the
different wall thicknesses of the hollow club head. The wall thicknesses generally vary
from around 1 mm on the top of the club, the crown, to 2.3 mm in the face of the club
head. All parts, except one, were modeled using the plastic kinematic hardening material
in LS-DYNA. This type of material declaration was used to model the high-strength
titanium used in most golf club drivers. The final part, which was the hosel where the
club head meets the shaft, was modeled as a rigid part. The club head was attached to the
shaft in the model by constraining two common nodes and, this material choice was
selected solely for joining the club head to the shaft. In the club head to golf ball impact,

the hosel does not experience much stress or deformation.

Table 7. Club head components.

 

 

Thickness
PART # Component of Clubhead Material (mm)
1 Face Titanium 2.41
2 Crown Titanium 1 .52
3 Sole Titanium 2.54
4 Backplate Titanium 1 .52
5 Hosel Steel 2.29
6 Heel Titanium 1 .52
7 Toe Titanium 1 .52
8 Face/Sole Interface Titanium 2.79

 

3.2.3 The “Arms”, “Wrists”, and “Hands”
The arms in the swing simulation were represented by a rather section of solid 8-
node elements. The elements were all modeled rigidly since in an actual golf robot these

elements are made of substantial steel parts. Their true purpose in the model was to

34

provide a foundation for the applied moments used to initiate and drive the motion of the
arms in the mechanical swing of the “Iron Byron”. The block representing the arms is
fixed at one end and a moment is applied to the center of mass.

Like the arms, the wrists were rigidly modeled using 8-node solid elements and
again their purpose was to provide a foundation for the applied moments in the model.
This block of elements was constrained to revolve around a node at the end of the arms.
The loading applied to this block of material helped represent the motion of the wrists in
the simulation.

The grip was also rigidly modeled and consisted of 4-node shell elements. On the
physical Iron Byrons, the grip is a collet mechanism that holds the grip of the club. The
shell elements contained in the grip surround the first 10 inches of the shaft in the model.
The grip played the role of the hands. It provided a location for the wrist-rolling moment
in the golf swing. A load was applied to the grip of the shaft, which allowed the golf club
to roll over towards the ball at impact, thus simulating the hands of the mechanical “Iron
Byron”.

3.2.4 The Swing: Loading, Joint Definitions, Contact Definitions

The loading in the golf swing consisted of 3 applied time varying moments
sufficient to independently drive the arm, wrist release, and wrist roll motions. As hinted
at above, one was applied to the rigid part acting as the arm, one to the part acting as the
wrists (wrist release), and one to the grip which acted as the hands (wrist roll). In the
finite element model, magnitudes of the moments were defined using 3 separate load
curves. These load curves had to be carefully iterated and tested to ensure that the swing

in the simulation resembled the swing of the mechanical “Iron Byron”. They were also

35

used to ensure that theiclub rotated correctly towards the ball making a relatively square
impact. The overall club head speed measured at impact in the swing simulation was
approximately 120 mph.

The joint definitions in the model also played a role in fine-tuning the swing by
allowing for a choice of the wrist release time as happens in the real golf swing. The
arms were connected to the wrists by creating a revolute joint in LS-DYNA. Therefore,
the part acting as the wrists can revolve around a certain point near the end of the arms.
Similarly, the end of the grip was also constrained to revolve around a certain point at the
wrists. At the beginning of the simulation the arms, wrists, and grip were all in a locked
position and the revolute joints were locked until certain times in the swing. This allows
the arms to start the overall downswing. As the block representing the arms rotates due
to the applied moment the revolute joint at the wrists unlocks and the club starts
revolving around the arms and the wrists. At a final time in the swing, the joint between
the wrists and grip unlock and the shaft rotates towards the ball just before impact.

The contact definition between the ball and club head became active in the model
when the clubface meets the ball. Both the static and dynamic coefficient of friction
values used were 0.3. To save computation time, the ball and club head were switched
from deformable to rigid bodies at the beginning of the simulation. Upon ball — club
head contact, the club head and ball switch back to deformable for the duration of impact.
Once the ball leaves the club head, both were switched. This rigid - deformable
switching dramatically cuts down on the computation time required for each simulation

of the swing.

36

3.2.5 The ball

The golf ball in the model consisted of 4 parts, all of which were meshed using 8-
node solid elements. The outer layer was modeled as a linear elastic material (moderate
stiffness ionomer), while the remaining inner layers were all modeled using the Mooney-
Rivlin-rubber material model in LS-DYNA. Essentially, this was a two-piece ball model
having single cover and core since the common nodes were merged at time of meshing.
The extra parts can be used to model more complex ball designs having cover, mantle,

and core(s).

37

3.3 HEEDS Optimization of Shaft Stiffness

3.3.1 Background / Aim

The optimization strategy presented in this part of the study hopes to establish a
method that not only evaluates the performance of various shaft designs but also aims to
discover unconventional shaft designs that may have been overlooked.

As shown in Chapter 2, the measure of a golf shaft’s “bending profile” or “E1
curve” is becoming a benchmark in the golf industry for the evaluation and comparison
of golf shaft designs. Some shafts are actually being designed to meet certain “bending
profiles”. The motivation for the optimization in this part of the study was to find
optimum “bending profiles” through evaluating a given full-swing finite element model.
The sections that follow explain the implementation of the optimization, its overall
objective function, as well as the results.

3.3.2 . Baseline Design

In order to prepare the model for the optimization, a traditional shaft design
needed to serve as a starting point. This starting point, or baseline design, did not need to
satisfy the constraints and objective. However, using a reasonable starting point helped
to ensure that the overall swing simulation was consistent and accurate. Details of the
problem’s constraints and objective will be explained in later sections.

The figures below show the baseline design at different points in the simulation.
The different points shown in the figures below are also used to compare the behavior of

the optimal designs found in the optimization.

38

 

 

Figure 12. Simulation at 0.15 seconds.

39

 

Figure 14. Simulation at 0.37 seconds.

40

._. .9

'-2Mo-.._,
.-. . .‘t
d.-

livfi.‘f)r
4“‘ Q
if a ..

.u. .- ‘
~‘-.-....'( . D .. 4
. _‘ ..
5.

Figure 15. Simulation at 0.41 seconds.

- f u‘. .. r,
. , l- .‘
r '1 ,- r:-

. m - o J '
misfiwmafl

 

Figure 16. Simulation at 0.42 seconds.

41

 

Q19 .7 *."‘\"

 

v \ 0
L. 8
Figure 18. Simulation at 0.44 seconds.

The baseline shaft was modeled after a standard steel golf shaft (True Temper

Dynamic Gold 8300). The shaft was cut into two inch segments and the radius and

thickness were measured. The moment of inertia could then be calculated for each

42

section of shaft. The figure below shows a breakdown of the stiffness profile used in the
baseline design. The range of stiffness in this shaft is typical of most steel and graphite

shafts in use today.

Baseline Shaft Design

El stiffness (N-mA2)

10.00 ~

 

 

o 9; 12113213233an;on
Length from butt (in)

Figure 19. Stiffness profile of baseline shaft design.

3.3.3 Preparations for Optimization
3 The original model of the swing simulation needed to be modified in various

ways to prepare for the optimization and to keep each evaluation as consistent as
possible. One of the initial problems involved keeping the impact of the golf ball on the
center of the clubface while the stiffness of the shaft was altered from run to run during
the optimization. Changes to the stiffness of the shaft would result in minor changes to
the location of the club head and location of the resulting impact with the golf ball.

The solution implemented was to separate the LS-DYNA run into two distinct
parts. The first part of the simulation was run without the golf ball model. The run was
then terminated once a predetermined node on the center of the clubface reached a

prescribed location. This location then became the location of impact with the ball. At

43

this termination point, a small Perl script was called that read the position of the
predetermined node on the clubface and then translated all nodes of the golf ball
accordingly to ensure that the ball made contact with the center of the clubface at impact.
This contributed to the consistency of each swing evaluation.

Having moved the ball to ensure it impacted the center of the clubface, the second
part of the simulation was run, which included the impact with the golf ball, and the
resulting launch and spin characteristics were then evaluated. A diagram detailing the
sequence of the simulation is shown below. The next section gives a detailed explanation

of the optimization problem.

 

 

1" Run

Golfball not included

Run terminates once clubface reaches pre-
determined lo cation

Location to become point of impact in 2"‘1 run

 

 

V

 

 

Translation of Golf Ball

Perl script translates nodes of golf ball to
coincide with the center of the clubface at the
termination location

 

 

 

2"ll Run

Golfball included
Run terminates shortly after impact
Rigid body data output

 

Figure 20. Sequence of swing simulation.

3.3.4 Optimization Problem Definition

To begin defining the optimization problem, the overall objective must first be
explained. As the stiffness and geometric properties of the shaft are changed in each run,
the resulting trajectory and spin characteristics of the golf ball are examined. The overall

objective is to minimize the golf ball’s backspin, while satisfying other important

constraints on trajectory, direction, and velocity magnitude.

45

 

 

 

As mentioned in a previous section, the swing speed in the simulation was near
120 mph. This swing speed is typical of most low-handicap or professional golfers. At
these high swing speeds too much backspin becomes an undesirable characteristic. High
backspin causes the golf ball to “balloon” or have too much lift, and can cause
considerable losses in driving distance. Therefore the overall objective of the
optimization was to find shaft designs that minimize backspin while subject to other

constraints on trajectory, direction, and velocity magnitude.

 

 

Figure 21. Variables influencing backspin.

The backspin and trajectory characteristics of a golf ball are functions of a limited
number of variables in the golf swing. The figure above depicts the three variables that
directly influence backspin; the angle of attack al, the clubface loft angle a;, and the club
head velocity V. These variables are also all functions of the dynamics of the golf shaft
during the downswing. In the downswing the shaft is in a bent back position and is

“kicked” forward towards the ball at impact. The degree of this action during the

46

downswing is determined by the bending stiffness properties of the shaft. The goal of the
optimization was to investigate the feasibility of unconventional shaft designs that meet
the problem objective and that satisfy the basic launch and spin constraints.

3.3.5 Model Variables, Responses, and Constraints

To reiterate, the goal during the optimization is to vary the bending stiffness
properties or E1 profile of the golf shaft. This was accomplished by keeping the Young’s
modulus, E, in the material definition constant while varying the moments of inertia Ixx
and Iyy in various elements in the LS-DYNA model. The moments of inertia [xx and Iyy
are functions of the shafts geometric dimensions, which are the outer radius R0 and
thickness I. These were the primary variables used in the optimization. The assumed
outer radius and thickness of each element were varied which in turn varied the moment
of inertia and the overall bending stiffness or E1 profile of the shaft.

In Section 3.1.1 the details of modeling the shaft in LS-DYN A were discussed
and must be noted here. Provided in this section is Table 6, which shows those elements
and sections varied during the optimization study. The total length of the shaft was 44
inches. However, portions of the butt and tip ends of the shaft, 12 inches and 3.8 inches
respectively, were kept at the standard baseline values. This was a result of the ends of
the shaft needing to meet certain dimensions in order to comply with the Rules of Golf
and to also allow for the universal assembly of grips and club heads.

The entire mid-section, comprising 28.2 inches, was varied in the optimization
study. Over this length the assumed outer radius and thickness for each element were
varied. These values were then used to compute the relative entries needed as inputs to

the LS-DYNA simulation. They were the moments of inertia bar and Iyy, the polar

47

moment of inertia J, and the cross-sectional area A. These were all functions of the
assumed outer radius and thickness of each element.

An important note to make is in regards to the torsional stiffness (similar to “El”
but termed “GJ”) of the golf shaft. The polar moment of inertia J is indirectly varied
during the optimization, which causes the torsional stiffness to vary as well. As stated
earlier, the focus of the optimization is on the bending properties of the shaft, however
the changes in the torsional stiffness must be addressed and acknowledged. It was
assumed that changes in torsional stiffness during the optimization would most
dramatically affect the dispersion (directional) angle of the ball flight. Therefore a rather
liberal constraint was put on the dispersion angle as will be explained later in this chapter.
The torsional stiffness of the shaft was not a major design consideration, however future
studies may want to include the modeling of composite golf shafts where this property
can be independently studied and optimized.

The difficulty in constraining the outer radius and thickness along the shaft will
now be discussed. In each design evaluation, the radii and thickness along the length of
the shaft must remain continuous and must also converge upon the butt and tip radii and
thickness as well.

The solution found was to use a combination of polynomials to define the radii
and thicknesses along the length of the shaft. This would provide continuity and also
allow for the values to converge to their tip and butt end values.

In order to implement this in the optimization, a small MATLAB program was
created to make these intermediate calculations. The inputs to the program were a series

of coordinate points that were used to define the polynomials and therefore the shape and

48

thickness of the shaft. Once the outer radii and thicknesses of the shaft were defined the
program’s outputs were then the moments of inertia Ixx and Iyy, which were then used as
inputs to the LS-DYNA swing simulation. The figures and equations below offer a
summary of this procedure. Example polynomials are given along with a table showing

the allowable range of the coordinate points in the Optimization.

Example Polynomial for Radius (n=2,3.4)

l I I T I

 

 

 

 

 

 

 

0,5 t. _
— Blended Polynomial
0.4 — Rm" -
--———— Rmin
75‘
g 0.3 .-‘_-__,_ -1 —~-—'— fi~~.\ 7
:6 \
s \

0.2 ~ \_

0.1r- ' —

 

 

 

0 1 I l l
0 5 10 15 20 25

Segment of shaft optimized (in)

Figure 22. An example polynomial for the outer radius of the shaft. The maximum
and minimum allowable radius is also shown.

49

Example Polynomial for Thickness (n=2,3)
0.04 . . f r .

—- Blended Polynomial (n=2,3)
" TmaX ' q
Tmin

 

 

0.035

T

 

 

 

 

0.03%.... ...... W

0.025 ~ 4

0.02 P / .4

Thickness (in)

0.015 . /’ -

OOOOOOOOOOOOOOOOOOOOOOOOOOOO

0.01 » _

0.005 » ~

 

 

 

0 5 10 15 20 25
Segment of shaft optimized (in)

Figure 23. An example polynomial for the thickness of the shaft. The maximum
and minimum allowable thickness is also shown.

The radius was defined using a combination of 2nd, 3rd, and 4th order polynomials
(n=2,3,4 in title of figure). The coefficients for each respective polynomial were found
and then the polynomials were summed and averaged. The generic polynomials below

help to explain how this was accomplished.

A=ax2+bx+c (3.3-1)
B = fx3 + gx2 + hx+ j (3.3-2)
C=kx4+mx3+nx2+px+q (3.3-3)

In order to solve for the coefficients of each polynomial, given the two end conditions,
each required 1, 2, and 3 coordinate points respectively. The thickness on the other hand
used a combination of only 2“d and 3rd order polynomials. Therefore the thickness

required 1 and 2 coordinate points for its solution. In all, defining the radius and

50

thickness along the length of the shaft required a total of 9 x and y coordinate points,
therefore 18 overall variables.

A couple of notes must be made regarding the choice of polynomials and the
choice of maximum and minimum radius and thickness. The use of higher ordered
polynomials was not considered due to the number of variables already being evaluated.
More design variables would significantly increase the design space and therefore the
number of runs in the optimization. The maximum and minimum constraints on the
radius and thickness may also seem conservative but were chosen to meet certain

feasibility and manufacturability requirements.

Table 8. Coordinates needed to define polynomials with ranges used in
optimization. The X coordinates range over the length of shaft being optimized (0-
28.2 in) while the Y coordinates range over the radius and thickness values
respectively (all units in inches).

 

X1 Y1 X2 Y2 X3 Y3

Radius

n=2 8 - 20 0.1 - 0.4 -- -- -- --

n=3 6 - 12 0.1 - 0.4 16 - 22 0.1 - 0.4 -- --

n=4 2 -8 0.1 -0.4 12- 16 0.1 -0.4 20-24 0.1 -0.4
Thickness

n=2 10 - 18 0.01 - 0.03 -- -- —- --

n=3 6 - 12 0.01 - 0.03 16 - 22 0.01 - 0.03 -- --

The formulas for calculation of the moments of inertia Ixx and lyy, polar moment
of inertia J, and area A of the shaft are given. The inner radius is simply found from

subtracting the thickness from the outer radius.

[xx = Iyy = at-(R: - ref) (3.34)

J = an: - ref) (3.35)

51

A = ”(R3 — Rf) (3.3-6)
The following diagram summarizes the entire sequence of each evaluation in the

optimization. This was similar to the previous diagram but also includes the MATLAB

implementation.

52

 

 

l

 

Poflomials in Matlab

0 Input variables are coordinates which define
polynomials for outer radius and thickness of shaft

0 Output are the moment of inertia’s Ixx, Iyy, the polar
moment of inertia .l, and cross-sectional areaA

needed in LS-DYNA simulation

 

 

 

l

 

1" LS-DYNA Run

0 Golf ball not included

Run terminates once clubface reaches predetermined
location

0 Location to become point of impact in 2"‘1 run

 

 

 

l

 

Translation of Golf Ball

0 Perl script translates nodes of golf ball to coincide
with the center of the clubface at the termination

 

 

 

location

 

2"“ LS-DYNA Run

Golf ball included

Run terminates shortly after impact

Rigid body data output

Evaluation of the ball flight and spin characteristics

 

 

 

 

 

Figure 24. Sequence of one evaluation in optimization.

53

The constraints and overall objective function of the optimization are now

discussed. Each run in the optimization was evaluated by the resulting golf ball flight

and spin characteristics after impact. Table 9 is a summary of the various constraints and

the objective function.

Table 9. Constraints and objective of optimization.

 

 

Response Type Target/Norm
X-vel Prerequisite

Y-vel Prerequisite

Z-vel Prerequisite

X-rot Prerequisite

Y-rot Prerequisite

Z-rot Prerequisite

Launch Angle Constraint >7 deg, <15 deg
Dispersion Angle Constraint <15 deg
Velocity Magnitude Constraint >140 mph
Y-rot Constraint <200 rad/sec
Backspin Objective Minimize

 

The formulas for calculating the launch and dispersion angles, velocity

magnitude, and backspin are given below. The x-axis is the main ball flight direction, the

z-axis the off-axis direction, and the y-axis the upward direction.

 

 

 

Launch Angle = arctan( MI (33.7)
vel
Dispersion Angle = arctan£ Zvel) (3.3-8)
Xvel
Velocity Magnitude = JXvel 2 + Yvel2 + Zvel 2 (3.3-9)
Backspin = Zrot(cos(DispAng )) -— Xrot(sin(DispAng )) (3.3-10)

54

3.3.6 Results and Evaluation of Optimal Designs

Once the problem was fully developed and defined in HEEDS, the optimization
run was executed. 140 different designs were evaluated during the optimization
sequence. These designs were then carefully examined.

In the initial stages of evaluating the results, those designs whose performance
numbers were similar in magnitude were kept while disregarding the rest. This left a
total of 22 shaft designs to be considered. The next step involved plotting their “EI
profiles” and eliminating many of the designs that were similar to each other. The “El
profiles” for the most unique designs will now be discussed and compared with the
baseline design.

Figure 25 shows the profiles of the first group of optimal designs. The number
assigned corresponds to the evaluation number in the optimization. Also shown is the
baseline design. As can be seen in Table 10 the overall “best” design found by HEEDS is
number 110. Designs that were close in performance but contained slight variations in

their profile are also shown in this group of designs.

55

Optimal Designs - Group 1

8.00E+01 -
7.008401 < .

60051-01 .

//

5.00E+01 <

4.00501 ~
7:) . ‘ ‘ #-M: 5;(7-.‘
3.00501 . ._ -—'- ""T‘h .

‘ fitmrrrfj‘cqu. --- 4 ~ __ N
2.00am . '

El stiffness (Nm"2)

1.00E+01 <

 

0.00E+00 . . . . .
0.00 5.00 10.00 15.00 20.00 25.00 30.00

Section of shaft varied in optimization (in)

 

 

[+110 +24 135 + 127 +Baseline Design]

 

Figure 25. Optimal Designs - Group 1. The profile is that of the 28.2 inch span
varied in the optimization.

Figure 26 contains the second group of optimal designs. These designs were
separated from the others due to their rather unique El profiles. They remain relatively
close to the “best” design in performance, but have rather unconventional profiles. Two
of the designs are similar in that their stiffness declines sharply and then contains a
distinctive increase in stiffness near the tip end of the shaft. These are very unique
designs but their feasibility must be questioned. The sharp decline in stiffness near the

butt end could cause the shaft to snap or break during the downswing.

56

Optimal Designs - Group 2

9006+01 1
8006+01 <
7.006001 <
6.wE+01 -
5.WE+01 -
4115401 <

Safe-01 -

 

El Stiffness (ch2)

2.1136401 .

Lflwl ‘

 

 

0.005000 - , . . . .
0.00 5.00 10.00 15.00 20.00 25.00 30 00

Section of shaft varied in Optimization (in)

 

|+es +61 65 +easoiine 011st

 

Figure 26. Optimal Designs - Group 2. The profile is that of the 28.2 inch span
varied in optimization.

Table 10 gives a summary of the major flight and spin characteristics for both
groups of designs. The velocities were all similar in magnitude and the directional angles
were well within the constrained limits in the optimization. However a note must be
made in regards to the launch angles and backspin rates. Clearly not all of the designs
included had necessarily the lowest backspin. This was due to a trade off made when
constraining the original optimization problem. The designs with the lowest backspin
generally had the lowest launch angles as well. This was somewhat undesirable. A
desirable launch angle is around 9-11 degrees, which is well known in the industry as
being the optimal launch angle for maximum distance. However it was felt that such a
narrow launch angle constraint would over constrain the optimization and therefore

overlook design considerations.

57

Table 10. Ball flight and spin characteristics of optimal designs.

 

Velocity Magnitude Directional Angle Launch Angle Backspin

 

Evaluation # Performance # (in/sec) (deg) (deg) (rpm)
1 10 0.85 3157 -6.48 7.99 2551
24 0.89 3171 -6.57 , 8.07 2674
135 0.92 3186 -2.62 8.21 2759
127 1.00 3197 1.28 8.54 2986
89 1.12 3216 -2.73 8.99 3362
61 1.17 3191 2.52 9.43 3516
65 2.63 3105 3.87 7.92 2296

 

Various designs above were re-run in LS-DYNA in order to compare their swings
with the baseline design and to make observations in regards to their feasibility. Each
design was chosen carefully.

When compared to the baseline design, no observable differences existed for the
first group of optimal designs, therefore making them good design considerations.
However the figure below shows design 65 as it was swung. The feasibility of the design
must be questioned due to the high degree of deformation occurring in the shaft. No

checking of stress levels was done in the current simulations.

saloon snow
The same A
it???
[23’
,5” if
13515,"

' w. . I“ _ «if

1-

\“

Figure 27. Design 65 exhibiting a high degree of deformation in downswing.

58

4 Conclusions and Discussion

4.1 Shaft Modeling

A majority of the experimental tests and methods used in the golf industry to help
characterize the performance of golf shafts can be modeled and studied analytically. This
was accomplished by an initial knowledge Of the shaft’s “bending profile” or “E1 curve”.
The “zonal frequency” tests as well as the various “kickpoint” tests were modeled using
simple beam finite elements.

The calculated zonal frequencies were over-predicted at the higher zones (shorter
cantilever lengths) but matched well at the lower zones (longer cantilever lengths). This
showed promise in eliminating the need for experimental measurements, but indicated
more investigation was needed into modeling the compliance of the frequency
measurement device.

The main conclusion drawn from modeling the “kickpoint” tests was that none of
these tests or methods lends any real insight into this somewhat elusive property of golf
shafts. These tests show no real distinction or relevance in their attempts of measuring
this property.

4.2 Optimization

An opportunity was taken to study and Optimize the golf shaft using a full swing

finite element model with golf ball impact. This allowed for a much broader study Of the

golf shaft and its influence in the golf swing. The actual dynamics Of the golf shaft and

its influence on ball flight and spin characteristics were studied and Optimized.

59

The results of the optimization revealed a rather unconventional group of shaft
designs as part of the designs considered. Consideration of stress limits on the shaft
material should be included in future Optimization runs. Prototyping Of these designs
should be pursued and testing established. As so many know in the golf industry, a new

unique shaft design could be marketed very easily.
4.3 Future Work

Based upon the results from the elementary shaft modeling and the successful
implementation Of the Optimization, there are many things to consider for future work.
Modeling Of Golf Shafts

0 Further investigation of modeling the Fujikura frequency machine and its

clamping mechanism

0 Researching new dynamic, rather than static, tests for characterization and

measurement of “kickpoint”
Optimgition

o Prototyping and testing Of the Optimal shaft designs

0 Shaft optimization for a range of swing speeds

o Re-examine loading and joint definitions to model more advanced robots

0 Examine experimental data of actual swing kinematics in hopes of modeling

an actual subjects golf swing

0 Incorporate modeling Of graphite composite golf shafts including ply-

orientation

0 Investigate using prescribed displacements and/or velocities rather than load

curves to gain more consistency in evaluating designs in the Optimization

60

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

Bibliography

Chou, A., and Roberts, O.C., “Golf Shaft Flex Point — An Analysis of
Measurement Techniques”, Science and Golf 11: Proceedings of the World
Scientific Congress of Golf, 1994, E&FN Spon, London, pp 278-283

Mase, T., “Correcting for local circumferential deflections when measuring golf
shaft flexural stiffness”, The Engineering of Sport 5 - Conference Proceedings of,
2004, M. Hubbard and RD. Mehta and J .M. Pallis, Chicago, pp 515-521

Horwood G.P., “Golf Shafts - a technical perspective”, Science and Golf 11:
Proceedings of the World Scientific Congress of Golf, 1994, E&FN Spon,
London, pp 247-258

Wallace, ES, and Hubbel, J E “The Effect of Golf Club Shaft Stiffness on Golf
Performance Variables -— Implications for Club-Fitting”, Materials and Science in
Sports, April 2001, pp 23-35

Butler, J .H., and Winfield, DC, “The dynamic performance Of the golf shaft
during the downswing”, Science and Golf 11: Proceedings of the World Scientific
Congress of Golf, 1994, E&FN Spon, London, pp 259-264

Brouillette M., “On Measuring the Flexural Rigidity Distribution of Golf Shafts”,
Science and Golf IV: Proceedings of the World Scientific Congress of Golf, 2002,
pp 387-401 '

Friedman, Z., and Kosmatka, J. B., “An improved two-node Timoshenko beam
finite element”, Computers & Structures, Volume 47, Issue 3, 1993, pp 473-481

Lee, M., and Kim, C., “Design Optimization Of graphite gold shafts based on
weight and dynamics of swing”, 5'“ International Conference on the Engineering
of Sport ( University of California, Davis), September 2004, pp 248-255

Reddy, J.N., An Introcuction to the Finite Element Method, 2ud Ed., New York:
McGraw Hill, 1993

Przemieniecki, J. 8., Theory of Matrix Structural Analysis, New York: McGraw
Hill, 1968

Cook R.D., Concepts and Applications of Finite Element Analysis, 2“d Ed., New
York: Wiley, 1981

Hughes, T. J .R., The Finite Element Method: Linear Static and Dynamic Finite
Element Analysis, Englewood Cliffs, N.J.: Prentice-Hall, 1987

61

[13] Han, S. M., Benaroya, H., and Wei, T., “Dynamics Of Transverse] y Vibrating
Beams Using Four Engineering Theories”, Journal of Sound and Vibration,
Volume 225, Issue 5, 2 September 1999, pp 935-988

[14] Winfield, DC, and Soriano, B.C, “Planar Motion of a Flexible Beam with a Tip

Mass Driven by Two Kinematic Rotational Degrees Of Freedom”, Journal of
Vibrations and Acoustics, Volume 120, January 1998, pp 206-213

62

 

 

 

 

 

 

 

 

   

 

   

   

 

   

 

 

 

 

 

 

 

 

   

 

 

 

 

 

     

   

 

 

   

         

   

 

 

 

 

 

     

 

 

 

 

   

   
 

 

  

l

'Tllilil

511’
l

      

6

2

 

‘7
i l
i

  

iii
01