Machine Tool Spindle Design PDF

Rochester Institute of Technology
RIT Scholar Works

Theses Thesis/Dissertation Collections
2000
Machine tool spindle design

Jamie Hoyt
Follow this and additional works at: http://scholarworks.rit.edu/theses
Recommended Citation
Hoyt, Jamie, "Machine tool spindle design" (2000). Thesis. Rochester Institute of Technology. Accessed from
This Thesis is brought to you for free and open access by the Thesis/Dissertation Collections at RIT Scholar Works. It has been accepted for inclusion
in Theses by an authorized administrator of RIT Scholar Works. For more information, please contact [email protected].
MAC~ETOOL
SPINDLE DESIGN
•Dy
Jamie A. Hoyt
A Thesis Submitted in
Partial Fuifiiiment of the
Requirements for i:he
MASTER OF SCiENCE
iN
MECHANICAL ENGiNEERING
Approved by:
Dr. Wayne W. Walter

Department ofMechanical Engineering Thesis Advisor)
Dr. Hany A. Ghoneim

Department ofMechanicai Engineering
Dr. Kevin B. Kochersberger

Department of Mechanical Engineering
Dr. Satish Kandiikar

Department ofMechanical Engineering
DEPARTMENT OF MECHANICAL ENGINEERING

ROCHESTER INSTITUTE OF TECHNOLOGY
MAY, 2000
Permission granted
Tide of Thesis: Machine Tool Spindle Design
i, iamie A. Hoyt, hereby grant permission to the Waiiace Library of the Rochester
Institute of Technoiogy to reproduce my thesis in whole or in part. Any reproduction will
not be for commercial use or profit.
Date: ~ •/17 /oc Signature of Author: _

Abstract
In modern machine tool applications the performance of a machine tool is judged
by its ability to produce work-pieces accurately and efficiently. The stiffness of the
machine tool spindle has a profound impact on the overall machine performance. The
work presented here provides a tool for machine tool spindle designers to develop
spindles that are sufficiently stiff to meet their needs. The analysis presented here is
divided into three main sections.
The first portion is a static analysis. The static analysis calculates the lateral
deflection of the spindle-bearing system. AMatlab program was developed that aiiows
the user to enter the spindle parameters into a batch file and obtain the plots of the
deformed shape of the spindle.
The next portion is a dynamic analysis of the spindle. This portion includes both
the modes of vibration and the forced response. The modal analysis treats the spindle as a
continuous Euier-Bernouiii beam. A numerical method for handling the steps in the shaft
and applied boundary conditions was developed that could be extended to many other
applications in rotor dynamics. AMatlab program was developed for the dynamic
analysis. This program provides a designer with plots of the mode shapes and forced
response given the spindle design parameters.
The final section is an optimization of the spindle design. Given constraints on
the location and stiffness of the support bearings, aMatlab program will return values for
these parameters resulting in the spindle configuration that presents the minimum
deflection at the spindle's gauge line.

Table of Contents
Page
Table of Contents 2
List of Tables 3
List of Figures 4
Nomenclature 6
1.0 Introduction 10
2.0 Static Analysis 15
2.1 Deflection of Elastic Shaft 17
2.2 Deflection of Bearings 26
2.3 Matlab Solution 27
3.0 Dynamic Analysis 34
3.1 Modal Analysis 35
3.2 Matlab Solution for Mode Shapes 61
3.3 Forced Response 66
3.4 Matlab Solution for Forced Response 70
4.0 Optimization Analysis "4
4. i Optimization Model 74
4.2 The Constrained Steepest Descent Algorithm 81
4.3 The Matlab Solution for Optimization 83
5.0 Conclusions 87
References 90
Appendix A : Batch Fiie Template 91
Appendix B: Matiab Programs for Static Analysis 93
Appendix C: Matiab Programs for Dynamic Analysis 107
Appendix D: Matiab Programs for Optimization Analysis 131

List of Tables
Page
3.1 Comparison of Resonant Frequencies (FEA vs. Analytical) 66
4.i Table of Design Variables 76
4.2 Optimum Values for Static Analysis 84
4.3 Optimum Values for Dynamic Analysis 84

List of Figures
Page
2.i Static Spindle Model 16
2.2a Elastic Deflection of Spindle Shaft 18
2.2b Deflection of Spindle Bearings 18
2.3 Model of Uniform Beam 19
2.4 Transformation of Beam Segment 22
2.5 Shear and Bending Moment Diagram for a Beam Segment 23
2.6 Batch File, Static Analysis 28
2.7 Malta Representation of Spindle Geometry 29
2.8a Deflection Contribution of Elastic Shaft 30
2.8b Deflection Contribution of Support Bearings 30
2.9 Total Deflection of Spindle 31
2.10 Comparison of Total Spindle Deflection (FEA vs. Malta) 33
3.1 Dynamic Spindle Model 36
3.2 Lateral Vibration of an Euler-Bernoulli Beam 37
3.3 Sample Euler-Bernoulli Beam 40
3.4 Mode Shapes for Sample Beam 43
3.5 Stepped Euler-Bernoulli Beam 45
3.6 Sample Stepped Euier-Bernouiii Beam 49
3.7a Boundary Conditions for Section i 52
3.7b Boundary Conditions for Section 2 52

3.7c Boundary Conditions for Section 3 53
3.7d Boundary Conditions for Section 4 53
3.8 Free Body Diagram, Joint 1 55
3.9 Free Body Diagram, Joint 2 57
3 iO
. Free Body Diagram, Joint 3 59
3.11 Batch File, Modal Analysis 62
3.12 Mode 1 Comparison 63
3.15 Batch File, Forced Response 7i
3.16 Comparison of Forced Response, FEA vs. Analytical 73
4.1 Optimization Model of Spindle 75
4. 1 Batch File, Optimization 85

Nomenclature
Symbol
A Cross sectional area of beam [in2]
A Gradient vector of inequality constraint [dependent on constraint]
ai Location of drive pulley [in]
a? Location of rear bearing [in]
a3 Location of front bearing [in]
&4 Location of gauge line [in]
bk Location of kth joint in spindle shaft [in]
c Gradient vector [in/in]
d Vector of design changes [unitless]
gi(x) Ith inequality constraint [dependent on contraint]
D Distance between support bearings [in]
E Young's modulus [psi]
f Quadratic subproblem [unitless]
FQ Cutting force [Ibf]
Fcub Unbalance force due to cutting tool [Ibf]
Fd Drive force [Ibf]
Fde Equivalent drive force [ibf]
Fdub Unbalance force due to unbalance of drive puiiey [ibf]
f(x) Cost function [in]
In Moment of inertia of nth beam segment [in4]

Kf Lateral stiffness of front bearing [Ibf/in]
Kfmax Maximum lateral stiffness of front bearing [ibf/in]
Kg Torsional stiffness of front bearing [in-ibfj
K; Generalized stiffness [ibf/in]
Kr Lateral stiffness of rear bearing [ibf/in]
Kfmax Maximum lateral stiffness of rear bearing [ibf/in]
M(x) Bending moment in spindle shaft [in-ibfj
MappHed Externally applied moments [in-ibfj
Mb Reaction moment at front bearing [in-lbfj
Mbe Equivalent reaction moment at front bearing [Ibf]
M; Generalized mass [lbf-s2/in]
Mk Moment induced at kth joint in spindle shaft [in-lbfj
Mt Moment at guage line due to cutting force and tool length [in-ibfj
OH Cantilever, distance between front bearing and gauge line [in]
q,; Generalized coordinate [in]
Qi Generalized force [Ibf]
R Penalty parameter [unitless]
Rf Reaction force at front bearing [ibf]
Rfe Equivalent reaction force at front bearing [Ibf]
Rr Reaction force at rear bearing [ibf]
Rre Equivalent reaction force at rear bearing [ibf]
Sj Siack variable for ith constraint [unitless]

T Kinetic energy [in-ibfj
tj Step size [dependent on design variable]
[Ty] Transformation matrix between ith and jth beam segments [unitless]
ti Length of cutting tool [in]
u Strain energy [in-lbfj
U Potential energy [in-ibfj
Uj Lagrange multiplier [unitless]
V Maximum constraint violation [dependent on constraint]
V(x) Shear force in spindle shaft [ibf]
Vk Shear force induced at kth joint in spindle shaft [ibf]
x Axial position along shaft [in]
x Vector of design variables
x0 Fraction of moment exerted

by front bearing [unitless]
yb Elastic deflection of spindle shaft [in]
y* Elastic deflection of spindle shaft [in. ]
6f Deflection at front bearing [in]
5q; Virtual displacement [in]
6r Deflection at rear bearing [in]
8wj Virtual work [in-lbfj
Ei Convergence criteria [unitless]
82 Maximum allowable constraint violation [unitless]
<j>i Ith normal mode

[unitless]
<p Descent function [in]
p Mass density [Ibf-sTin4]
to Circular frequency [rad/sec]

1,0 Introduction:
Great demands are placed on the capabilities of today's modern machine tools to
produce parts that are dimensionaiiy correct with increasing accuracy and throughput.
Some of the machine tool components that impact the accuracy and throughput of the
machine are the drive systems, way systems, control and feedback systems, and finally
the machine tool spindle. The machine tool spindle is the element of the machine that
either supports the work-piece or the cutting tool. In addition to being a support
structure, the spindle also rotates at high rates of speed to provide relative motion
between the work-piece and the cutting tool. Therefore the spindle has a direct impact on
both the throughput (material removal rate), and the accuracy of the finished part.
According to Lewinschai (1985), the most common requirements of a machine
tool spindle are:
High running accuracy
High speed capability
Great stiffness
Low and even running temperature
Minimum need of maintenance
Often in machine tool spindles these parameters will conflict with each other. In order to
achieve a higher speed capability the designer must trade off spindle stiffness for speed or
visa-
versa. The spindle designer must carefully weigh the requirements of the user to
determine the best possible balance of these parameters.
10
The goal of this research is to provide a tool for a spindle designer to aid in the
evaluation of the spindle stiffness. High running accuracy, high operating speed
capability, low and even running temperature, and minimum need of maintenance are
typically functions of the bearing's geometry, manufacturing, lubrication, and method of
mounting. If the spindle designer is able to quantify the stiffness requirements for the
bearing he can then work with the bearing manufacturer to select the proper bearings for
the application.
Al-Shareef et ai. (1990) developed a quasi-static method of analyzing machine
tool spindles. Their analysis takes the amplitude of the dynamic forces and applies them
to a static model of the spindie-bearing system. For the static analysis the deflection
contribution of the spindle shaft and the deflection contribution of the spindle support
bearings are superimposed to obtain the total deflection of the system.
The static analysis of the spindle shaft assumes that a stepped flexible shaft is
pinned in the location of the support bearings. The analysis of this flexible shaft consists
of a transformation from a stepped shaft to a uniform shaft. This transformation yielded
additional shear and

bending moments at each of the joints in the shaft. The resulting
uniform shaft was analyzed classical mechanics.

using
The deflection contribution of the spindle support bearings assumes a rigid shaft
supported by linear springs. The reaction forces yielded the deflection at each of the
springs. Essentially, the deflection contribution of the bearings is a straight line fit
between the resulting deformed positions of the springs.
11
In addition to the static analysis an optimization of the deflection at the end of the
spindie was presented. The optimization analysis consisted

primarily of varying the
spindle design parameters and looking at the effect on the resulting deflection at the
spindie gauge line. Plots were presented iiiustrating the effect of the variation of these
parameters. The following conclusions were drawn from these plots.
In the design of a spindie there exist an optimum ratio of the bearing spacing
to the overhang of the spindie. As the flexurai stiffness increases and the ratio
of front to rear bearing stiffness decreases the optimum bearing-overhang
decreases.
A dimenskmless flexurai stiffness (Kf(OH)3/EI) of greater than 5 results in
minimum deflection at the cutting tool. The deflection at the end, or gauge
line, of the spindie is very sensitive to the flexurai stiffness for magnitudes
less than 5.
Having more than 3 steps in the shaft is desirable for obtaining minimal
deflection values.
The magnitude, position, and direction of the driving force greatly effects the
deflection at the gauge iine. For each scenario there exists an optimum
location of the drive pulley.
In Lewinchai (1983) a similar study on the variation of spindie design parameters
was presented. Plots were generated that illustrated the effect of the bearing spacing-
overhang ratio on the spindie stiffness for support bearings of varying stiffness. From
these plots it could be concluded that for very stiff support bearings the optimum
spacing
12
between the bearings becomes shorter. It couid also be concluded that if the spindie has a
long overhang the stiffness of the bearings has a lesser impact on the stiffness of the
spindle.
Other work in the optimum design of machine tool spindles was also done in
Montusiewicz et ai. (1997). In this work a model of a machine tooi spindie supported by
hydrostatic bearings was presented. Tne study consisted of applying a four-stage
muiticriterion optimization strategy to a static modei of a spindie. The objective of the
analysis was to reduce the radial and axial deflection of a spindie, the total mass of the
spindie, the totai power loss of the bearings, and finaiiy the size of the bearings. The
analysis divides the spindie system into four subsystems. Each of these systems are
optimized locally, and finaiiy integrated to provide a giobal optimization. The outcome
of this analysis was a computer aided optimum design package. This package aiiows
spindie designers to interactiveiy design an optimum spindie, inputting required design
variables throughout the optimization process.
A qualitative dynamic analysis of a machine tooi spindie was presented in Al-
Shareef et ai. (i99I). Traditionally in the dynamic analysis of machine tooi spindies the
first mode is thought to be responsible for poor cutting quality. The purpose of this work
was to assess this assumption. There was concern that this would not be the case since
the range of operating frequencies for a given spindle often excite the higher modes. The
first four modes for an example spindie were solved for analytically and compared to
experimentai results. The modal analysis presented ignores damping and rotational
affects. The authors site an experimentai study that proved there to be little difference
between the non-rotationai natural frequencies and the rotational criticai speeds. By
looking at the individual mode shapes they found that the first mode contributed the most
to the deflection at the tooi to work-piece interface. Aii other modes in the operating
frequency range exhibited nodal characteristics at this interface. Since the excitation
force would be exerted here they concluded that the first mode would indeed be most
accountable for poor cutting quaiity. However they also noted that at the higher modes
there was significant deflection at the location of the support bearings. This couid result
in the degradation of these bearing and an eventual loss of spindle stiffness.
Some other works, pertaining more generally to the field of rotor dynamics, were
also researched. Two of these works deal primarily with the extension of the conventional
transformation matrix (CTM) technique. In the work done by Curti et ai. (1993) an
expression for an 8 x 8 dynamic stiffness matrix of a rotating Timoshenko beam is
derived and reiated to the conventional 4x4 dynamic stiffness matrix. This provides for
the inclusion of anisotropic supports.
In work done by Murphy (1993) a polynomial transfer matrix was deveioped to
replace the conventional transfer matrix for modal and forced response analyses. The
advantage of the polynomial transfer matrix is an increase in computational speed of 3.5
to 100 times over the conventional transfer matrix. Example problems were analyzed
using both the CTM and PTM methods as weii as a finite element analysis. The results
for ail three cases were identical and the speed of the PTM method was considerably
faster.
14
2.0 Static Anaiysis:
The static anaiysis calculates the lateral deflection of the spindie. Figure 2.1
illustrates the model under scrutiny. The following assumptions were necessary to
perform the anaiysis:
1 . The spindie shaft is assumed to be an Euier-Bernouiii Beam.
2. The spindie is subjected to a cutting force, a drive force, and the reaction forces at the
bearings. The drive force must be applied behind the rear bearing.
3. The torsional and axial deflections of the spindie shaft are neglected.
4. The centeriine of the spindie shaft is exactly iniine with the centeriine of the bearing
bores. There is no contribution to the lateral deflection due to manufacturing
misalignment.
5. The spindie housing and the cutting tooi are both assumed to have an infinite
stiffness.
6. It is assumed that the spindie is supported by only two bearings. This is common for
most machine tool spindles. Manufacturability precludes the use of more than two
bearings in most spindies.
7. The contribution of transverse shear deformation to the overall lateral deflection is
assumed to be negligible. It was observed in a study conducted by Al-Shareef and
Brandon, that the contribution of shear deformation is dependent on the ratio between
the length of the spindie and the spindie nose overhang. The shear deflection for
short spindies with small overhangs contributes more to the overall
15
Figure 2. i Spindie Model
16
deflection than longer, more siender spindies. A variety of spindies were analyzed in this
study and a maximum contribution of 12 per cent was found (Ai-Shareef et ai., 1990)
Superposition was employed to calculate the iaterai deflection of the spindie. The
elastic deformation of the spindie shaft, ys and the deflection of the spindie bearings, yb
were superimposed to calculate the overall deflection of the spindie (see figs. 2.2a and
2.2b). Equation 2.1 gives the overall deflection of the spindie.
v. =ys+yb (2.1)
2.1 Deformation of Elastic Shaft:
For the elastic contribution of the spindie shaft Ai-Shareef and Brandon propose a
method to transform the stepped spindie shaft to a uniform shaft (Al-Shareef et. ai, 1990).
This approach will be employed in this analysis. When the shaft is transformed there is a
moment, M* and shear force, \\ induced at each step in the shaft (fig. 2.3). In addition
the applied forces and reactions must be transformed into equivalent forces applied to
beam segments with larger bending moments of inertia. These equivalent forces are
"e"
noted using the subcript (i.e. Fj -> F<je).
The deflection of the uniform beam can be easily analyzed using conventional
beam theory and singularity functions. The singularity functions will be represented by
expressions in < >. If the value of the expression within these brackets is less than zero
<2-4>2
the function becomes zero (i.e. =

0). If the value of the expression is greater than
<4-2>"
zero, the function simply becomes the expression within the brackets (i.e. =
<4-2f).
The shear force, V(x) of the uniform beam can be found to be:
17
a1 H|
Figure 2.2a Elastic Deflection of Spindie Shaft
'////
Figure 2.2b Deflection of Spindie Bearings
18
1 X] Oi XI OJ \1 b4
O* XI
m b5
to XI DO XI TOW. XI tJIl-1 X] 311 X
*
11 1 1~1 1 1
VI
i_
V2i V3, V4i V5i V6, Vn-2, Vi
^M2^M3^M4^M^Mr^
1
I
Fde
XSxVx
Rre
X>xV\.
Rfe
Mbe
"J
a2 y a3 a4+ti
Figure 2.3 Model of Uniform Beam
19
F(x-u1)D
+Fc(x-a,-t
Vix) =
(2.2)
-Rfe(x-a>r+iv*(*-b*r
The moment of the beam, M(x) becomes:
M(x) =
]V(x)dx+MappI,eJ (2.3)
Fde{x-a1)1
M(x) = + Fc{x-a4 -ttf -R^x-aJ

W
+
Vt{x-bkf +Mk(x-bky
*=i
The slope of the beam, 8(x) becomes:

i
0(x) = --JM(x)dx
(2.5)
^{x-aiY+^(x-at-tiy-^(x-a^-^(x-a,y
0(x) =
(2.6)
El
+M{x-<*2) +1i
Integrating the slope of the beam yields the elastic deflection, ys(x):
{t-f(x-a$'6
o
+!f{x-a4 -tlY -%<x-a,V
o
'
y,(*) =
-^(x-atf+tîx-bj+t^x-btf (2.7)
El. o i=i 6 i=1 2
+
Mu , ,2
+qlx +
-f-(x-a3) q2
The integration constants, qi & q2 can be found by applying the following boundary
conditions:
ys(x =
a2) = v
20
Solving for the integration constants yields:
o2 =
~%a2 -a;f
6
-
*=i
^-(a2 -bkf
o
-^{a2 -bkf -qxa2
k=l 2
(2.8)
1
<fe
fc3 ~{2 -i>3)-TL{3

6
(2.9)
(2-3) t=i 6
+
I S
f^L-b.)2-(a,-b.)2)
2
The derivation for the moments and shear forces induced, and the equivalent
applied forces when the stepped shaft is transformed into a uniform shaft will now be
presented. The derivation begins by looking at the internal shear and bending moments
for an arbitrary segment in the stepped beam (Fig 2.4). An illustration of the shear and
bending moment diagrams is also offered (Fig. 2.5).
From the shear and bending moment diagrams it was found that:
V{x) =
Vl=Vr (2.10)
and
M(x)=M,-V,x (2.11)
From Castighano's Second Theorem:
dU
7=y (2.12)
m
and
*>;r t
ou
6 (2.13)
dM
21
The strain energy, U for one-dimensional
bending is known to be:
'7777777777777777777777777777777'/. NJ/
Figure 2.4 Transformation of a Beam Segment
22
V
/lx
J Vr
Figure 2.5 Shear and Bending Moment Diagrams for a Beam Segment
23
2
r M
H^n* <214>
o
It should be noted that this expression for the strain energy does not include any
contribution due to transverse shear deformation. Substituting equation (2. 14) into
equations (2. 12) and (2.13), with eqns. (2. 10) and (2. 1 1) for the original beam segment
(prior to the transformation to the uniform shaft), yields the following y and 6:
*
EI\ 2 3 J
i f y p }
e =
\Mrl--^-\ (2.16)
^
2 }
El\
Similarly, when the analysis is repeated for the segment after its transformation the
8*
deflection and slope,
y*
and are found to be:
'
v*=-\MZ-J.2Li
EI* (2.17)
2 3 j
{
0* = -\\MJ^-\
r (2.18)
EI 2
[ J
The differences in y and 8 must be compensated for with the induced shear force and
bending moment.
(2.19)
[EI EI 2 3 )
}\
l
A^ ^J Umi-LtL-).
=
f2.20)
^
[/ / 2 J
J[
Therefore:
24
\\MJ>
i 1 1 \}Mf Vf\_\ 1
VJT] (2-21)
[EI 2 6
J 2 6
ET)\ [ET}{ j
md (2.22)
[EI EI )[ 2
J [EI 2
}[ J
Multiplying eqn. (2.22) by -1/2 and adding eqn. (2.21) yields:
1
V =/ 1_ (2.23)
/ r
Substituting eqn. (2.23) into eqn. (2.21) yields:
1 J.
MM=r ?w (2.24)
i r
This anaiysis can be repeated to find an induced shear and bending moment at
each step in the shaft. The induced force and moment now become applied forces to a
beam segment with a moment of inertia of I This analysis can be extended to show that
aii applied forces must be scaled by a factor of In/I. Where In is the moment of inertia of
the uniform beam (largest moment of inertia in the stepped shaft), and I is the moment of
inertia of the segment that the force is appiied to. Therefore the induced forces and
moments become:
1 i
vk=i 7~
(2-25)
r
"i i
Mk=In -~]M
(2.26)
where:
) ) >*
Vr =Rr{x- +Rf{x-
fd{x fc{x (2.27)
ttf
- - - -
a, a4
-
a2 a3
25
Mr =RXbk-a2Xbk-a2y+Rf(bk-a3)(bk-a3y
-fc{K-<**-ti)(bk-a4-tl)-mb(bk-ay
The cutting, driving, and reaction forces from the stepped spindie shaft must also
be scaled to provide equivalent applied forces on the uniform shaft. The scaling of these
forces yieids:
E*=i~Fd (2.29)
Ifi
K=-rK (2 30)
*Rr
*Rf
Mbe=-^Mb (2.32)
2.2 Deflection of Bearings:
The deflection contribution of the spindle bearings was calculated by assuming
that the spindie is a rigid shaft supported by two flexible bearings (Figure 2.2b). The
cutting force Fc, and the driving force Fd were used to soive for the reactions at the
bearing. The two reaction forces were used to calculate 6r and 5f, the deflections at
the two bearings. The deflection contribution of the bearings is a straight line through
8r and 8f.
(Sr-Sflx-a2)+Sr(a3-a2)
yt kû)
-
t \
{<*3-<*z)
26
Where
8 -
fcfo +a-a*)-mh-fA<>* ~a2))
(234)
(a3 -a2)Kr
Sf
fAas -a2) + mb-fc(p4 +tl-a3)-(fc + fd\a3_a2)
(a3-a2)Kf (2.35)
mb=fMA+d-<**K (2 36)
2.3 Matlab Solution:
A program was developed using Matlab to automate the static anaiysis of the
spindie shaft. The user must simply enter the geometry, loads, and support parameters
file"
into a spreadsheet called a iSbatch A copy of the batch fiie template is presented in
Appendix A. The Matlab programming code used to automate the static anaiysis can be
found in Appendix B.
An exampie of the anaiysis for a simpie spindie is presented here. Figure 2.6
illustrates the batch fiie for the static anaiysis. Upon the completion of the batch fiie the
program wiii read the fiie and report a geometric representation of the spindle. The plot
illustrates the geometry of the shaft as well as the locations of the bearings, cutting force,
and drive force (see figure 2.7). This feedback allows the user to easily check for
mistakes in the batch fiie. With ail the information correct the program calculates and
reports plots of the deflection contribution of the elastic shaft (figure 2.8a) and the
deflection contribution of the bearings (figure 2.8b). Finaiiy the program reports a piot of
the totai deformation of the spindle (see figure 2.9).
27
Batch File:
Geometry:
Number of Sectjons(#):
Section Length OuterDiameter InnerDiameter Area Moment of Inertia

(#) fin-) (ml (inl fin "21 fin "41
1 3 2.25 2 0.83449 0.47265783
2 3 2.375 2.125 0.88357 0.560861726
3 3 2.5 2.25 0.93266 0.659419991
4 3 2.625 2.375 0.98175 0 76890787
5 3 2.75 2.5 103084 0.889900605
6 3 2.875 2625 1.07992 1022973438
7 0 0
8 0 0
9 0 0
10 0 0
11 0 0
12 0 0
13 0 0
14 0 0
15 0 0
Bearings:
Lat Stiffness of Rear Bearing QbAri): 100000

Lat Stiffness of Front:Bearing Qb/n): 500000
Fraction of mom. on Front Bearing; 0.1

Location of Rear Bearing (in.) 75
Location of Front Bearing (in.) 13.5
Pulley:
Location of Pulley (in.)
Static Belt Tension (lb):

~x-nK$M-v<>:..-'.;;:
Tool:
'i^mlî^M^^^M: "zM^ffi&ai&m,
Static Cutting Force fib): j 300
W^^^^^^t^^^- I 330
t^^^-fÛî^f'^^-.-^j'
*&>*-?
^33:33-
-','''.-
^mmSi^timWî *#>^flB*5SS4S3i
Length of Tool (in) 1 2
Speed:
Material Prooerties:
Modulus of Elasticity (psi):
Figure 2.6 Batch File for Matlab Solution
28
input dimensions
1 1 1 1 T 1 1 1 1
OD
-
ID
? D Pulley
7 -
-e- -v-
Rear Bearing
+ -*-
Front Bearing
~
Tool
6 -
-5-
4 -
3 -
...
r
* 1
,_,
1 1 1 1 1 1 1
8 10 12 14 16 18
x,(in)
Figure 2.7 Matlab Representation of Geometry
29
Deflection Contribution of Elastic Shaft
x 10
Figure 2.8a Deflection Contribution of Elastic Shaft
Oeflection Contribution of Bearings

x 10
1
-i
Figure 2.8b Deflection Contribution of Support Bearings
30
Combined Spindle Deflection
x 10
Figure 2.9 Total Deflection of Spindle
31
In order to confirm the results offered by the program, a finite element analysis of
the sample spindle was using Ansys. The

one-
performed spindle was modeled using
dimensional linearly elastic beam elements. The bearings were modeled using linear
spring elements. The cutting force was transformed into a force moment couple and
applied at the end of the spindle shaft in order to account for the tool length. Figure 2. 1 0
compares the deflections of the shaft using both methods. It is clear from the plot that
there is an excellent correlation between the two analyses.
32
Static Deflection Comparison
1.00E-02
8.00E-03
? FEA
=- O.OOE+00 +
>- Matlab
*
X(in)
Figure 2. 10 Comparison of Total Spindle Deflection (FEA vs. Matlab)
33
5.0 Dynamic Analysis:
The dynamic anaiysis for the spindie shaft consists of two portions. The first part
of the analysis is the modai analysis. The beam is treated as a continuous system for this
portion of the anaiysis. The second part of the anaiysis solves for the deflection of the
spindie by means of modal superposition. The following assumptions were made in
order to perform the anaiysis:
1 . The spindie shaft is assumed to be an Euler-Bernoulli Beam.
2. The spindie is subjected to a cutting force (FcSin(rj>ct)), a drive force (FdSin(Qdt)),
unbalance forces (FcubSin(cot) & (FdubSin(ot)), and the reaction forces at the bearings.
The drive force must be applied behind the rear bearing. The cutting force and drive
force are assumed to be harmonic.
3. The masses of the puiiey and cutting tooi are assumed to be concentrated. The mass
of the puiiey is assumed to be concentrated at the centeriine of the puiiey. The mass
of the tooi is assumed to be concentrated at the end of the spindie shaft. This point is
often referred to as the gauge iine.
4. There is no unbalance excitation introduced by the spindle shaft.
5. The rotational affects of the spindie shaft are neglected.
6. The torsionai and axiai deflections of the spindle shaft are neglected.
7. The centeriine of the spindle shaft is exactly iniine with the centeriine of the bearing
bores. There is no contribution to the lateral deflection due to manufacturing
misalignment.
34
8. The spindie housing and the cutting tooi are both assumed to have an infinite
stiffness.
9. It is assumed that the spindie is supported by only two bearings. This is common for
most machine tooi spindles. Manufacturabiiity typically precludes the use of more
than two bearings in most spindies.
10. The contribution of transverse shear deformation to the overall lateral deflection is
assumed to be negligible.
1 i. Damping is neglected in the dynamic anaiysis.
The modei scrutinized in the dynamic anaiysis is very simiiar to the modei used in the
static analysis. One major difference is the use of a torsional spring to represent the
torsional stiffness of the front support bearings. In addition the masses of the puiiey and
cutting tooi are included. See figure 3. 1 for the dynamic modei under scrutiny.
3.1 Modal Anaiysis:
The foundation for the modal anaiysis is the derivation of the wave equation for
the iaterai vibration of a continuous Euler-Bernoulli beam. Figure 3.2 represents the free
body diagram of an differential element of an E-B beam. Applying Newton's second iaw
to the beam element it can be shown that:
~^
=
-pA~^ (31)
dx at
and
V^ (3.2)
OX
35
bn
bk
b2
3
b1 ^
Mt
FdSin(wdt)
Md
"X 12 Tool i _~.
FcSin(wct)
Kr
FcubSin(wt)
FdubSin<wt)
/?77777 mrrw fmrm
a3
a4 TL
Figure 3.1 Dynamic Spindie Modei
36
y
">
Figure 3.2 Differential Element of an E-B Beam
37
It can also be shown from Strengths of Materials that:
M =
EI^4 (33)
Substituting eq. 3.3 into 3.2 yieids:
V =
EI^4 ox-
(3.4)
Finaliy, substituting eq. 3.4 into 3.1 and rearranging yields:
d"v hJ d"v
-+^r^r
dx*
= 0 0.5)
dr pA
The following harmonic solution to eq. 3.5 was assumed:
y(x,t) = y(x)smo)t (3.6)
Substituting the assumed solution (eq. 3.6), into the differential equation (3.5) yieids the
following forth-order differential equation:
dx4 P*y = 0 (3.7)
where:
P=B*EL
p (3.8)
EI
It can be shown that the general solution to the preceding forth-order differential equation
is:
y(x) = A cosh fix + B sinh fix + C cos fix + D sin fix (3 .9)
38
Equation 3.9 represents the wave equation for an E-B beam. The mode shapes for a
beam can be found by substituting values for p that correspond to the resonant
frequencies. The constants A,B,C, and D can be solved for by applying the boundary
condition for the beam.
A systematic method involving numerical methods was developed to soive for the
resonant frequencies and their corresponding mode shapes. This method is not exclusive
to the spindie probiem at hand. It can be extended to the lateral vibration of many Euier-
Bernoulii problems. Listed below are the steps to this method:
1 . Establish the boundary conditions for the system.
2. Collect the system of equations into matrix form.
3 .
Using Gaussian Elimination numerically reduce the matrix.
4. Using the Bisection Method or a comparable root finding method soive for the
resonant frequency, p.
5 . Back substitute to find the constants A,B,C and D for the beam segment.
Figure 3.3 represents a simple beam used to illustrate this approach. Tne beam under
scrutiny here is a uniform E-B beam fixed at both ends. The first step is to find the
boundary conditions. Since the beam is fixed-fixed, the displacement and rotation at x
=
0,i are both equal to zero. Expressed mathematicaiiy:
K0) = 0 (3.10)
7(0) =
0 (3.11)
y(i)=v ^3.iz)
7(0-0 (3.13)
39
x=u
Figure 3.3 Sample Euier-Bernouiii Beam

Substituting Eq. 3.9 into Eqs. 3.10-3.13 yieids:
V(U) = A +C = 0 V-iV
y [V) =
B + jj =
v
y(l)
-
a cosn(p/) + d smn(/) +c
cos(/?/) + u sm( fit) =
v (j.io;
y (/)
=
A sinh( />/) + B cosh(/?/) -
C sin( fii) + D cos(/w) =

0 (3 17)
The next step is to collect this system of four equations into matrix form. This yieids Eq.
3.18:
10 10 \A 0
0 10 1 B 0
< ? = >
(3.18)
cosh(/3?) sinh(yS0 cos(/#) sin(/?/) 0
sinh(>9/) cosh(y9/) -sin(yS7) cos(/37) > 0
Guass- Jordan elimination. The

Step (3) reduces the matrix in eqn. 3. 18 using reduced
system is illustrated in eqn. 3.19:
1 0 1 0 rA"
V
0 1 0 1 B 0
0 0 1
sin(/ff)-sinh(/ff)
. = <
o} (3.19)
cosOGf/)-cosh03T)
0 0 0 cos0G0cosh(y3/)-l K.
DJ L
0J
The reduced system can be used to soive for the resonant frequencies, pj
(cos(fil)cosh(fil) -l)D = 0 (3.20)
If D was equai to zero, then AB, and C would aiso equal zero. This would not be a
meaningful result. Therefore it can be concluded that:
(cos(j0/)coshOSZ)-l) = O (3.21)
This is where the root finding method suggested in step (4) comes into piace. The
roots of eq. 3.21 iead to the resonant frequencies of the system. Solving for the roots
yieids:
fiti,fi2l,fi3l =
4.7,7.8.1 1.0
After solving for the roots the final step is to back substitute to obtain the
constants A3,C and D. Begin the substitution by assuming that D=i Working .
backward from D it can be shown that the remaining constants are:
sinh(ffi)-sin(/7/)
cos(p/)-cosnij&/)
B = -l
.
_
sin(/ff)-sinh(/?/)
-
cos(/) cosh(/?/)
Substituting these constants into eq. 3.9 yieids the mode shape for the sample beam. The
equation for the mode shapes becomes:
sin(/y)-sinh(/?y/)
cos(B 7l) -
coshifi} ,1)
J
j=U,.. (3.22)
cosifijl) cosh(yff;7)
Figure 3.4 illustrates the first three mode shapes for the sample beam.
The method described here can be applied to find the mode shapes for ail uniform
E-B beam problems. However if the beam is stepped, as is the case with the spindie
shaft, there needs to be a set of boundary conditions for each beam segment. This leads
42
First Three modes for Sample Seam
mode 1
- - - -
mode 2
mode 3
V //
-0.5
\
\
w>
-1.5
x (norm alized)
Figure 3.4 Mode Shapes for Sample Beam
43
to a very iarge system of equations. Aside from the problem of having a very large
system, the number of steps wouid change for different spindies. This wouid make
automation very difficult. A transformation matrix was developed to handle the steps in
the shaft. The transformation matrix relates the constants on one side of a step to the
constants on the other side of the step. This makes the number of equations in the system
independent of the number of steps in the shaft.
The deveiopment of this transformation matrix begins by looking at an arbitrary
step in an Euier-Bernouiii beam (see figure 3.5). In order for continuity to exist the
deflection, slope, moment, and shear force at the joint must be the same for both beam
segments.
=
yi\i) y2Kl) v>-<;
JV(/) =
Jy(/) (3-24)
{EIl^Ml = {Ell^Ml (3.25)

ax" ax"
mm={EI),qM
'
0.26)
ax ax
Substituting eq. 3.8 yieids:
yjj smh(yDj/) + C, cos(y&,/)

=
ax cosn(pj/) +
-+-
uK sm(p,/)
(3.27)
A2 cosh(/y) + B2 sinh(/y) + C2 cos(J32l) + D2 sin(fi2l)
fit(A sinK^,/) + 5, cosh(/9/) C, sin(#/) + Dx cosifij))

- =
(3.28)
fi2 (A2 sinh(/y) + B2 cosh(j32l) C2 sm(fi2l) + D2 cos(B2l))
-
fiÂ, cosh(#/) + Br sinh(#/) -C, cos^fij) Dx sin(#0) =

-
(3.29)
fi2{A2 cosh(fi2l)+B2 sinh(/?2/)-C2 cos(fi2l)-D2 sin(/?2/))
44
3
\
\ \
El El
Figure 3.5 Step in E-B Beam
45
/V(4 sinh(#/) +
fi, cosh{#/) +C, sin(#/)
-
>, cos(#/)) =
(3,30)
/?23(^2 smh(&/) +2 cosh(p2l) + C2 sin(fi2l) ~D2 cos0?2/))
The system of four equations and eight unknowns can be coiiected into matrix form
0091(0,;) smh(/y) cos(/y) sm(AO -cosh(/M) -sdnfa(p,i) -cos(p,/)
smh(#/) codtfyW) -saAfi,}) ooe<J,l) -smh(J?,/) -cosb^,/) -oos(JJ)

, 2
(&\0i (Ethfii (EJ\J: (EI\P?

- - -r- -r-
- -
snuî) siring) o^l^)

-
srrinî) coshil,!) -cosf,^,ij ooA^fi^j smtn^j 3

(7).tf. (EI\fii {EI\B{ (EI).^
'a 0
A 0
Q 0
A 0
.=.
0
(331)
5, 0
n
c2
iA,
0
Using Gauss-Jordan elimination foiiowed by back substitution a relationship can be
found between Ai-Di and A2-D2. Two ratios, Ri and n, were defined to simplify the
reiationsmp.
(^h
^ =
(3.32)
(0,
1=^
(3.33)
46
A, =
tl/?1
+1i(cosh0g1/)cos(^2/)-/?1 sinh(#/)siiih(2/)K
\rR2
+\]
+ -* '{cosKÔsin^/?,/)-^, sinh(#/)cosh(y92/)}82
(3.34)
_ K _
^{coshOgjOcosOg,/) +
R} sinKAOshX^/)^
Cli?1
-
~1i{eo8hOg1Osin(^2/)-^l sinKAOcoÔK
(3.35)
+**-*
2
'{, co&(0tl)w(B2l) + sinhÔcosÔJCj
+
E1^1
2
J{sinh()g/)sirt(/y) -
i?, cosh(BJ)cos(B2l)p2
tli?1
C, =
'fa sin(#/) sinh(/y) -
cos(y9I/)cosh(y92/)}^2
+ M_ziJ^?] si^B1l)cosh(fi2l)-&nhifi1l)cos(fi2r)}B2
(3.36)
+
^ +1^
{cos(#/)cos<2/) + fl, sin(#/) sm(/y)}C2
+-l^^^s(^1/)sin(^/)-i?Jcos{^1/)sin<y92/)}D2
47
Vr2
-ii
i-J J
-
{#, cos09,/)cosh(/?2/) + sin(#/) sinh(2/)}52

L 2 1 (337)
-
~^Rl2+%l coÔshX^/) sin</y)sin</y))C2
+
*li?1
'{ft, coÔcosOff,/) + sin(#/) sin(/?2/)}D2
The coefficients from eqns. 3.34-3.37 can be collected into a transformation matrix [Tj,
such that:
'
Ax 'K
B2
< =
{Tj >
(3.38)
IA, D,
The use of the transformation matrix can be iiiustrated by expanding the sample
beam problem to inciude steps in the beam (see figure 3.6). Applying the boundary
conditions wouid result in the following system of equations:
'K
By 0
10 10 0 0 0 0 0
1C)
0 10 10 0 0 0 0
< = <
(3.39)
0 0 0 0 eosh0ff3/3) sinh(^3/3) cos(B3l3) sm(fi3l3)
}A] 0
0 0 0 0 0
sinh(&/3) cosh(/?3/3) -sin(j03/3) cos(jff3/3)
\B>] 0
c3
0
A; * J
If transformation matrices were not used, the only way to solve the system of equations
wouid be to relate A>-Di to A3-D3 by including the continuity equations. This wouid
increase the size of the system to 12 equations and 12 unknowns. It wouid also make the
48
Figure 3.6 Sampie Stepped Euier-Bernouili Beam
49
size or tne system aepenaent on tne nunroer or steps in tne snan. 1 nis in-turn wouia
make automation more difficult. If the transformation matrices were used the system of
equations wouid be reduced to 4 equations and 4 unknowns, regardiess of the number of
steps in the shaft. The first two equations in the system become:
10 10 B3
0 10 1
K\T] > =
< (3.40)
c3
The last two equations wiii be the same as represented in eqn. 3.39. Once the system of
equations is deveioped steps 3-5 of the pre-described method can be used to solve for the
resonant frequencies and their corresponding mode shapes.
The five-step process and transformation matrix can now be combined and
appiied to find the frequencies and modes shapes of the spindle depicted in Figure 3.1.
In order to encompass aii of the

externally applied boundary conditions the beam must be
divided into four sections. Figures 3.7a-3.7d depict the four subdivisions. The first
section is between the rear free end and the drive puiiey. The second section is between
the puiiey and the rear support bearing. The third section is between the rear and front
support bearings. The forth and final section is between the front support bearing and the
cutting tooi. There wiii be four constants for each of the four sections for a total of
sixteen constants.
Beginning with the free end of section one, the shear force and bending moment
at x
=
0 are both equal to zero.
Expressed mathematicaiiy:
50
EI^- 0
V,(0) = =
(3.41)
ax
and
EI^-
M,(0) = = 0 (3.42)
Substituting eqn. 3.9 into equations 3.10 and 3. 1 1 and setting x equal to zero yieids:
5,-D, =0
(3.43)
and
?~ s\
At the junction between sections 1 and 2 there are four boundary conditions. The
first three conditions involve the deflection, slope and bending moment at the joint
between sections i and 2. Since there are no externally applied moments, and the
structure is continuous, the deflection, slope, and bending moment at the joint must be
equai for both sections. Therefore:
Mai)=yiUh) (345)
JV(,) =
JV(a.) (346>
Eiyx \a1)
=
niy2 \a,) (j.<w;
Substituting equation 3.9 into equations 3.45-3.47 yieids:
A, cosh(/5b, ) + Bi sinh( Sal ) + C5 cos(pat ) + Dt sin( fial )

(3.48)
-
A2 coshOfibj ) B2 sinh( fiax ) C2 cos(fiax ) D2 sin( fiax )

- - -
= 0
At sinh(/x7j ) -+-
B. cosh(y(5bi )
C, sin(/3ar. ) + D, cosOSa, )
-
(3.49)
-
,42 sinh( Bax)-B2 coshfjSa, ) + C2 sin( fiax)-D2 cos(J3ax )

=
0
51
VI
Md /
/
V2
a1
4 3
Figure 3.7a Boundary Conditions for Section 1
V2(a1)
ms*\<
Md
a1
///////
Figure 3.7b Boundary Conditions for Section 2

a^
V3(a3)
V3(a2) \1'
V4(a3)
a3
Kr X Kf.Ktt
/////////
Figure 3.7c Boundary Conditions for Section 3
Mt
V4(a3)
a3
Kf.Ktt
/////////
Figure 3.7d Boundary Conditions for Section 4
53
Ax cosh(/xjt) + .81 sinh(/x*5)-C5 cosOSa.)-D, sin^)
-
A2 coshfjfifatj ) B2 sinh(y8afj ) + C2 cosO^ ) + D2 sm(Bax )

-
= 0
The forth boundary condition at this joint is affected by the mass of the puiiey. The mass
of the puiiey introduces an external shear force. Figure 3.8 iiiustrates the free body
D'
diagram at the joint. The shear force introduced by the mass is equal to the Aiembert
force associated with the puiiey mass.
Therefore:
Vm =
mpy =
-mdm2y2(ax ) (3.51)
For equilibrium at the joint:
rr y \ -w t \ rr /** J*>\
/
yâi)-y2(ql)
=
vm (i.oz)
Ely^(ax)-Ely2^(ax) =
-mda}2y2(ax) (3.53)
Substituting equation 3.9 into equation 3.53:
Ai sinh(/w, ) + Bx cosh(pax ) + C5 sinf/Sz, ) D, cosGSa, ) -
-4[sinKÎ)-^^cosh(^0]-57[cosh(/fa1)-^^sinh(/b1)] -
(3.54)
p EI fi'EI
+C2[sin(Î)--f
^^cos(fiax
fi3EI
cos^^)] '
+ Atcos^)--!
Afoostfb,
2
x )
sm(p\)]
-
^
fi3EI
= 0
The first three boundary conditions for the joint between the second and third
section are the same as the boundary conditions between the first and second joint.
Therefore:
A2 cosh( fi(a2 a, )) + B2 swh(fi(a2 ax )) + C2 cos(fi(a2 a, ))

-
- -
+ D2 sin( fi(a2 ax )) A3 cosh(fi{a2 ax )) B3 sinh( fi(a2 a, )) (3.55)

- -
- - -
C3 cos(/?(a2 -
ax ))
-
D3 sin( fi(a2
-
ax ))
= 0
54
V1(a1)
V2(a1) /
Vm%
Figure 3.8 Free Body Diagram of Joint 1
35
a2 %\wi\p\a2 -ai))+n2 cemj>[a2 ~ax))-K.2 $m(p(a2 -a, jj
+ D2 cos(fi(a2 -ax))-
A3 sinh( fi\a2
-
ax ))
-
B3 cosh(/5(a2
-
ax )) (3.56)
+
C3 sm(fi(a2 -
ax ))
-
D3 cos(B(a2
-
a.
)) = 0
A2 cosh(p'(a 2 -
a, )) + 52 sinh(p(a2
-
a, ))
-
C2 cos(p(a2 -
a, ))
-D2 sin(/?(a2 -a,))-/43 cosh(/?(a2 -a,))-53 sinh(p(a2 -ax)) (3.57)
+ C3 cos(fi(a2 -ax)) + D3 sin(p(a2 -
a, ))
=
0
For the forth boundary condition at this joint the shear force introduced by the rear
support bearing must be accounted. Figure 3.9 iiiustrates the free body diagram at the
joint. The shear force introduced by the bearing is proportionaito the shaft's
displacement at the joint.
r^-
rr s \ ?<* ^C\
vkr=Kry3ia2) (xrt)
rr / \ rr s \ rr s*y J?C\\
EIy2!,,(a2)-EIy3"!(a2) =
KTy2(a2) (3.60)
Substituting equation 3.9 into equation 3.60:
A2 sinh( fi(a2 B2 cosh( fi\a2 -ax)) + C2 sin(p(a2 ax ))

-
-ax))+
-D2 cos(B(a2 -ax))-A3[swh(fi(a2 -ax)) + -^-cosh(fi(a2 -a,))]

fi'EI
K* U
-53[cosh(/?(a2 -ax)) + -Smh(B(a2 -ax))]+C3[sm(B(a2 -ax))
fi'EI
+-^-cos(>9(a2 -aI))]+
D3[cos(B(a2 -ax))^-sin(fi(a2 -ax))]
= 0
fi EI fi EI
The first two boundary conditions for the joint between the third and fourth
section are comprised of the

continuity conditions (y3=y4 and y3'=y4?).
56
\ /*i / o\
vâ2)
V3(a2).
Vkr si/
rigure 3.9 Free Body Diagram of Joint 2
57
Therefore.
A2 cosh(/?(a2 -
a, )) + B2 sinh(p(a2
-
ax )) + C2 cos(/?(a2
-
a, ))
+ D2 sin(p(a2 cosh(/?(a2 (3.62)
-ax))-A3 -aj)-53 sinh(/>(a2 -a,))
-
C3 cos(/>(a2 -ax))-
D3 sin( fi(a2 -ax))

=
Q
A2 sinh(/>(a2 a, )) + B2 cosh(fi(a2 a, )) C2 sin(^(a2 a, ))

- - - -
+D2 cos(B(a2 -ax))-A3 sinh(/?(a2 -ax))-B3 cosh(/J(a2 -a,)) (3.63)

+
C3 sin(/?(a2
-
ax ))
-
D3 cos(/?(a2 -
a5 )) = 0
The third and forth boundary conditions are influenced by the bending moment and
shear force associated with the torsional and iaterai stiffness of the front support bearing.
Figure 3.10 illustrates the free body diagram at joint 3.
s~
rr / \ rr / \ rr f+ *~\
ia3 ) v j od;
=
^v3 ^3 )
-
r,iy4 a ,
v4 ia3 ;
and
/
i "* \ * / ^ /I jCiCi
^3l3;-Aj4(a3;=iWtr/ lJOOJ
EIy%"{a3) -EIy4'(a3) =
Kfy4(<*3) (3 67)
Substituting eq. 3.9 into eqs. 3.65 and 3.67 yieids.
A3 sinh(/?(a3 -a2))+B3 cosh(pi>3 -a2)) + C3 sin(^(a3 -a2))
-D3 cos(B(a3 -a2))-A4[smh(B(a3 -a2))+^-cosh(>9(a3 -a2))]
ATy .,,,.
in ^ (3-68)
54[cosh(^(a3 -a2)) + ^-sinh(y9(a3 -
a2 ))] + C4[sin(y?(a3 -a2))
+
^-cos(^(3
LI
-2))]+>4[cos(>9(^ '
-a,))-^-sin(>9(a, -a,))]
EI
= 0
fi fi
58
M3(a3)
V3(a3)
M4(a3) V4(a3)
\y Vkf
Mkft
Figure 3. 10 Free Body Diagram of Joint 3
59
ana
A3 cosh(p (a3 -a2)) + B3 sinh( fi (a3 a2 )) C3 cos(/?(a3 a2 ))

- - -
-A sin(/?(fl3 -a,))-4[cosh(/?(a3 -a2))+^-sinh<Aa3 -a,))]
^'^
-54[anh(#a, -fll)) + -^7Cosh(>9(fl7 -a^M+CJcos^^, -a,))
+ -^sinO?(a,
"
-a,))]+
Z)4[sra(jS(a, -a7))--^-cos(/?(a,
" -a,))]
= 0
/?'/ B"EI
The final two boundary conditions are related to the cutting end of the spindie.
The first of these conditions relates to the bending moment. Since the rotary inertia of
the cutting tooi is neglected the moment at the end of the spindie is equal to zero.
Therefore:
EI^-^- 0
Ml(0) =
dx"
=
(3.70)
A4 cosh(p(a4 -
a3 )) + B4 sinh(p(a4
-
a3 ))
-
C4 cos(fi(a4 -a3))-
D4 sin(/?(<z4
-
a3 ))
=
0
The last boundary condition involves the shear force at the end of the shaft. The shear
force is equai to the D'Aiembert force associated with the mass of the tooi.
Therfore:
V(a4 ) =
mty(a4 )
=
-mp 2y(a4 ) (3 .72)
Substituting equation 3.9:
60
A4[smh(P(a4
-
a3 )) + ^-cosh(y?(a4
El
-
a3 ))] + Bs[cosh(8(a4
-
a, ))
fi
2 2
+ ^-sinh(/?(a4
EI
-a3))] + C4[sin(y?(a4 -aJ) + %-cas(P(a4
fi'EI
-a,))] (3.73)
fi
m,(o"
.
r-*
-D4[cos(fi(a4-a3))^wi(P(at-a3))]
i =
0
fi'EI
In order to soive the set of simuitaneous equations, the set of sixteen equations and
sixteen unknowns were collected into a matrix.
3.2 Matiab Soiution for Mode Shapes:
A program was developed using Matlab to automate the modal anaiysis of the
spindie shaft. Data is coiiected and entered into a spreadsheet. This sheet acts as the
batch fiie for the modai analysis. Much like the static anaiysis, the user must enter the
geometry, mass information, and support parameters into the batch fiie. A copy of the
batch fiie template is presented in Appendix A The Matlab programming code used to
automate the modai anaiysis can be found in Appendix C.
An example of the anaiysis for a simpie spindie is presented here. Figure 311
iiiustrates the batch fiie for the

our"
modai anaiysis. The "grayed information does not
pertain to the modal analysis. Upon the completion of the batch file the program will
read the fiie and report a geometric representation of the information. With aii the
information correct the program calculates and reports the resonant frequencies for the
sample spindie. The sample spindie was also modeled using Ansys. A comparison
between the FEA and analytical results for the first three modes is presented in Figures
3. 12 through 3. 14.
61
Batch File:
Geometry:
Number of Sections(#): ZJ
Section Length Outer Diameter Inner Diameter Area Moment of Inertia
(#) On) On) On) 0n~2) 0n~4)
1 3 2.25 2 083449 0.47265783
2 3 2.375 2.125 088357 0.560861726
3 3 2.5 2.25 0.93266 0.659419991
4 3 2625 2.375 098175 0.76890787
5 3 2.75 2 5 1 03084 0.889900605
6 3 2.875 2.625 1 07992 1 .022973438
7 0 0
8 0 0
9 0 0
10 0 0
11 0 0
12 0 0
13 0 0
14 0 0
15 0 0
Bearings:
Lat. Stiffness Rear

Bearing (lb/In):
of 100000
Lat Stiffness Front
Bearing (lb/in):
of 500000
Tor^tiffnessofFrontBearing On- lb)
Location of Rear Bearing (in.) 7.5

Pulley:
Location of Pulley (in.) 4.5

Mass of (Ib-^gVinfr
Pulley
J
m
i 33.33
emmmmsb
Tool:
0 025879917
^ttR
MaRaseSc^:uJi^:."ftwe*>fe.); :
Length of Tool (in) | 2 I

Speed:
Material Properties:
Modulus of Elasticity (psi): 30000000
Density (lb^n~3): 0.289
Figure 3.11 Batch File for Sample Spindle
62
Mode 1 Comparison
(FEA vs. Closed Form)
-
-0.2
-0.4
X (in)
Figure 3.12 Mode 1 Comparison
63
Mode 2 Comparison
0.9
0.8 -
0.7
0.6
>-
0.5 -
0.4
0.3
CF Solution
0.2 FEA Solution
0.1
10 15 20
X (in)
64
Mode 3 Comparison
0.8 -
0.6 -
0.4 -
0.2 -
2<?
-
-0.2
-0.4
CF Solution
-0.6 FEA Solution
-0.8
X (in)
65
The resonant frequencies for the first three modes are compared in table 3.1. It is
clear from the table that the two methods correlate very closely for the two methods.
Table 3.1 Comparison of Resonant Frequencies (FEA vs. Analytical)
Mode FEA Analytical Difference

j
(Hz) (Hz) (%)
1 19.26 19.17 0.4672897
2 61.28 61.72 0.7180157
3 95.57 100.33 4.9806425
3.3 Forced Response:
The forced response of the spindle is calculated using a numeric modal
summation procedure. The development of the forced response begins with the equation
of motion for a beam, Dahleh et. aL (1989).
[/v"(x,of +m(x)y(x,t) =
f(x,t) (3.73)
The normal modes for the beam, <|>i(x), must satisfy the following equation:
(EI6)"-a>Mx)t'i=0 (3.74)
In addition to eqn. 3.74, since the normal modes are orthogonal they must also satisfy the
following equation:
i&rfjdx -
0 for * j (3.75)
0
The solution to the forced response can be represented in terms <j>,(x) as:
y(x,t) =
2>,(x)<7,(0 (3.76)
66
Where q,(t) is the generalized coordinate. The generalized coordinate can be realized
using the Lagrange Equation. Looking first at the kinetic energy yields:
T = \]y\x,t)m{x)dx
2
0
Substituting eqn. 3.76 for y(x,t) yields:
T =
^HHaiaj
>JMj*"(x)dx
2', ;
T =
\ZMa (377>
Where the generalized mass, Mj is defined as:
Ml=^2(x)m(x)dx (3.78)
o
The potential energy, U can be defined as:
U = ÊIy"2(x,t)dx
2 0
1 ' J 0
^4w (379)
2 .
Where the generalized stiffness, Kj equals
K,=\EI[<p"(x)]2ax (3.80)
0
If eqn. 3.74 is substituted into eqn. 3.79 it can be shown that:
67
tf =
(3.81)
A generalized force, Qi can be defined by looking at the work done by a virtual
displacement, 6qi.
<*,
=
J/(*,o2>^*
rearranging:
<^,=2>,a (3.82)
where:
ft=J/(x,fM(*>* (3.83)
From the Lagrange Equation:
f dT\ 8T dU
(3.84)
dt
Substituting for the kinetic energy, potential energy, and the generalized force yields the
following differential equation.
j/(*,0*,(x>fc
<ii+a>i2<li=T (3.85)
j<f,2(x)m(x)dx
where:
. (3.86)
68
The model of the spindle assumes four simple harmonic loads. The harmonic
loads include the drive force, cutting force, unbalance of the pulley and unbalance of the
cutting tool (see fig. 3.1). All four of the forces are assumed to be in phase with each
other and of the form:
f(x,t) =
F(x)sm(cot) (3.87)
Each of the forces are applied to a single point. Assuming the force is applied at x =
x,, it
can be described using the delta dirac function as:
f(x,t) =
Fsin(a>t)S(x-x0) (3.88)
By definition the delta dirac function is equal to zero for all x not equal to Xo. Further it
can be shown that:
]F(x)S(x-xo)dx =
F(x0) (3.89)
Substituting this relationship into eqn. 3.85 yields:
4t^=4WW
(3.90)
j<f>2(x)m(x)dx
Assuming the following solution to eqn. 3.90:
qi(t) =
q,sin(0Jt) (3.91)
yields:
ft
=
\co2
M"f (392)
-G)2)\<t>2mdx
69
The denominator of eqn. 3.92 must be broken down for the four sections of the spindle
and each of the segments (steps) in the shaft described in the modal analysis.
a, = tl^K (3.93)
W _ft)2E'wS J (Ak cosîx)+Bk sinh09,x)+Ct cos(0ix)+Dk
sin(#x))2
dx
n=\ k=\ o
For this analysis only the summation of the first four modes were utilized. After
the first four modes the difference between the resonant frequencies and the drive
frequencies become large and q, approaches zero. Therefore the steady state response
becomes:
Y =
Mx + <t>2<l2 +
&<73 + Ma (3 -94)
The deflections, Y were calculated for each of the four excitation forces and superposed
to yield the total forced response:
y,=y*+Ym+Y*+YA* (3-95>
3.4 Matlab Solution for Forced Response:
A program was developed using Matlab to automate the calculation of the forced
response for the spindle shaft. The magnitude and frequency of the excitation forces is
entered into a batch file. In addition to the load information the program reads the first
four modes calculated in the modal analysis program. The Matlab programming code
used to automate the forced response can be found in Appendix C.
An example of the analysis for a simple spindle is presented here. Figure 3.15
information
out"
illustrates the batch file used for this example problem. The "grayed
does not pertain to this analysis. It should be noted that the program will not function
70
Batch File:
Geometry:
Number of Sectjons(#): 3
Of) (in) On) fin) (in~2) fin-4)
1 3 a25 2 0.834486 047265783
2 3 2.375 2.125 0.883573 0.560861726
3 3 2.5 Z25 093266 0.659419991
4 3 2.625 2.375 0.981748 0.76890787
5 3 2.75 2.5 1030835 0.889900605
6 3 2.875 2 625 1079922 1.022973438
7 0 0
8 0 0
9 0 0
10 0 0
11 0 0
12 0 0
13 0 0
14 0 0
15 0 0
Bearings:
Lat Stiffness of Rear Bearing (ItMn): 10000C

Lat Stiffness of FrontBearing (lb/In.):
Tor Stiffness of Front Bearing (in- lb):
Location of RearBeanng(in)
Location of Front
Bearing (in.)
Pullev:
Location of Pulley (in.)
Mass of Pulley (\b-s~2Ari
Harmonic Drive Force (lb):

Drive Frequency (Hz):
Pulley Unbalance (rb-s~2):
Tool:
Mass of Tool Oj>^n):
Harmonic Cutting Force (lb): 30C

Cutting Frequency (Hz): 133.33
Tool Unbalance (lb-s"2): 0.004559453
Length of Tool (in)
Speed:
Spindle Shaft Speed (Hz): 16.666665671

Density (lb/fn'3): 0.289
Figure 3.15 Batch File for Sample Spindle Forced Response
71
properly if the modal analysis from section 3.2 is not completed first. The sample spindle
was also modeled using Ansys. A comparison between the FEA and analytical results for
the forced response is presented in Figure 3.16. The comparison between the FEA and
analytical responses shows a close correlation between the two methods. There is a 6.5%
difference in the deflection at the tool for the two methods.
72
Dynamic Response Comparison
0.01
0.008 ->-
0.006 +
0.004 4-
0.002 +
FEA Response
-Analytical Response
-0.002
-0.004
-0.006
-0.008
-0.01
X (in.)
Figure 3.16 Comparison of Forced Response FEA vs. Analytical
73
4.0 Optimization Anaiysis:
The optimization anaiysis consists of

minimizing the deflection of the spindie
shaft at the gauge line (see figure 4. 1). This anaiysis builds upon the static as weii as the
dynamic anaiysis. Optimal parameters are offered for both cases. The foiiowing
assumptions appiy to the optimization anaiysis:
i . The design variables for this anaiysis are the iaterai stiffness and the position of the
two bearings. All other parameters are assumed to be constant.
2. Each design iteration is approximated using the Tayior series expansion. This
approximation is required to define a quadratic programming subprobiem.
3. The optimization point may or may not be the global minimum. However the vaiues
assure a local minimum.
4.i Optimization Modei:
The development of the optimization problem rests in minimizing a cost function,
f(x), where x is the design variable vector. For the optimization of the machine tooi
spindie the cost function, f is defined as the deflection at the spindle's gauge line.
/(i) =
^(a4) (4i>
Given vaiues for the design parameters, a value for yt(a4) can be obtained numerically
using the Matlab routines developed in Chapters 2& 3. The design variabies, x are listed
in table 4. 1 . The remainder of the spindie design parameters are assumed to be fixed.
This is a fairly accurate assessment since for an existing spindie design the other
parameters wouid significantly influence the supporting components (i.e. gearbox and
spindie housing).
74
Figure 4.1 Optimization Modei of Spinaie
Tabie 4.1 Table of Design Variables
Design Variable Vector Parameter

xfl) a(2), postion of rear org
x(2) a(3), position of front brg
x(3) 1 Kf, lateral stiffness of front bra.
x(4) j Kr, lateral stiffness of rear brg
Generai constrained optimum design defines the following equality and inequality
constraints respectively:
g,.(x)<u v^.j;
For this optimization problem there exists no equality constraints. The following
equations define the inequality constraints.
xx>ax+D (4.4)
x2<a4-GH (4.5)
x3<KfimK V)
*4 <
Krnax (4.7)
To summarize these constraints, the first constraint (eqn. 4.4) stipulates that the location
of the rear bearing must be beyond the location of the puiiey by a distance, D. This is
"outboard"
required to ensure that the pulley is of the support bearings and there is
sufficient spacing to accommodate the width of the puiiey and the width of the bearing.
The second constraint (eqn. 4.5) requires that there exist a sufficient overhang to
accommodate features in the spindie shaft to accept and support the tooi. Tne third and
bearings'
forth constraints (eqns. 4.6-4.7) ensure that the stiffness vaiues are physically
obtainable. Without these constraints the optimization could potentiaiiy specify a bearing
with an infinite lateral stiffness.
Prior to developing the process used for this optimization anaiysis it is important
to first introduce the Lagrange function and the Lagrange Multiplier Theorem, Arora,
(1989). For general constrained optimization the form of the Lagrange equation is:
n m
1=1
Since there are no equality constraints in this anaiysis the Lagrange equation reduces to:
i=l
From eqn 4.9, f\x) is the cost function, m is the number of constraint equations, uis the
f1 i"1
lagrange mulitplier for the constraint equation, gi is the constraint equation, and Sj the
t"
slack variable for the constraint equation. The slack variable is a constant that
converts the inequality constraint to an equality constraint.
gj(x) + S2=0 (4.10)
If the design point, x is a locai minimum the Lagrange Muiitpiier Theorem stipulates the
following Kuhn-Tucker Conditions:
3L
=
0 forj =
lton (4.11)
dxj
dL
= 0 fori=ltom (4.12)
CL,
= 0 fori=ltom (4.13)
els.
77
f" f*5
If the constraint is inacti%'e Si is equal to zero. If the constraint is active Hi is equal to
zero. Therefore:
u1si
= 0 (4.14)
The Lagrange multipliers and the slack variables can be found by solving this system of
equations (eqns 4. 1 1-4. 14).
In order to apply numerical methods to soive for the design change an
approximate quadratic programming subprobiem (QP subprobiem) was defined. The QP
subprobiem can be obtained from a Taylor series expansion of the cost function. It has a
quadratic cost function and linear constraints. The probiem is defined as the
minimization of:
5(d*d/'
/ =
c*d'*'
+ 0. (4. 1 5)
where:
f(xw dw /<xw
/ = + ) -
) (4. 16)
subject to the following constraints:
A*dw<p (4.17)
/t-\ *i-
km
where
<r*";
is a vector of changes in the design variables for the design point, c is a
vector containing the gradient of the cost function f^x^*), and A is the gradient of the
inequality constraints.
In order to soive the QP probiem a search direction and a step size must be
determined. The constrained steepest descent method was used to soive for these two
entities. When no constraints exist the search direction is simpiy in the direction of the
78
negative of the gradient vector (d=-c). In the case of the spindie optimization constraints
exist and they must be included in the development of a search direction.
In order to accommodate the constraints in solving for a search direction a descent
function must be defined for the constrained probiem. A descent function must possess
two properties. First, it must be equai to the cost function at the optimum point. Next it
must allow for a unit step size near the optimum point. This is important because a unit
step size wiii yieid a high rate of convergence. The Pshenichny's descent function <p was
chosen since it obeys these two rules.
<pix; i <t 1 o;
/ tx; + a v ix,
=
In eqn 4. 18, R is the penalty parameter and V is the maximum constraint violation. The
user specifies the initiai vaiue of R. A subsequent value for R is calculated at the end
each iteration in the optimization process. In order to satisfy the necessary condition the
penalty parameter must be greater than or equai to the sum of Lagrange multipliers at the
k*
iteration.
R>rk (4.19)
For the m constraint equations, rk is defined as:
r,*=2>/ (42)
j=i
k*
for the
1th
is
Ujk
where the Lagrange multiplier constraint at the design point. The
Lagrange muitipiiers can be found by solving the system of equations

previously
mentioned (eqns. 4. 1 1-4. 14).
kfe
The maximum constraint violation at the iteration, \\ is defined as:
Vk=TMXp,gx,g2,....,gm} (421)
The next step in solving the optimization probiem is to define a step size
determination procedure. The decent function wiii yield the search direction, the step
size determination wiii dictate how far to adjust the design variables in that direction. For
this analysis an inexact line search method was used. For this method a sequence of triai
step sizes, tj was defined.
(1V / A r\r%\
t.=\-\ forj=0.1.2 (4.22)
Each iteration begins with the triai step size to=i. If a defined descent condition is not
satisfied the step size is cut in naif (ti=I/2). For a step size iteration, j and a descent
iteration, k the new design variable vector is defined as:
x<*+u)
=
xi*)
+
^dw
(4.23)
The acceptable step size wiii be the smallest integer j that satisfies the descent condition.
Qw^t-'A (424>
where *k+i j is the descent function defined in eqn, 4, 1 i evaluated at the trial step size.
d^7
The constant fik is found using the search direction,
fik=rYk)\ (4-25)
The constant y is specified by the user and has a value between 0 and 1 . The value of y
affects the allowable step size. Larger vaiues of y wiii result in smaller vaiues for the step
size. The end result is a slower rate of convergence. Alternatively very small vaiues for
80
y can iead to instabilities in the optimization process. Typicaiiy experimentation takes
place to find a suitable value for the engineering problem being solved.
iterations of search direction and step size are continued untii the method
converges on a iocai minimum for the cost function. Convergence is defined as the
design point were
|jdjj<^ (4.26)
where Ei is a specified small positive number.
4.2 Tne Constrained Steepest Descent Aigorithm:
A CSD aigorithm was used to optimize the design variables in a spindie shaft.
This section describes the steps to this aigorithm.
The first step to the CSD aigorithm is to set the counter, k equai to zero. At this
step initial vaiues for the design variables x, the penalty parameter R, the constant y, and
the convergence criteria si. An additional convergence criteria was aiso added to the
anaiysis. It was stipulated that the maximum constraint violation, Vk must not exceed a
predefined vaiue 82. This assures that design points with excessive constraint violations
are not allowed. A vaiue for this constant is aiso needed at this step. Since the goal of
this anaiysis is to optimize an existing spindie design the initial vaiues for the design
variables would simply be the parameters used in the existing design. The initial value
for the penalty parameter was defined as R=i . The constant y was defined as 0.5. Finaiiy
the vaiues for the convergence constants Si and 82 were both defined to be 0. 1 .
81
The next step is to calculate vaiues for the cost and constraint functions as weii as
their gradients. It is important to note here that the design variables were normaiized for
this anaiysis. Since the magnitudes be

of the variables
vary significantly it wouid
inappropriate to use there gradients in obtaining a search directions and step size. The
gradient of the cost function was calculated by applying the forward difference method to
the static and dynamic modei previously developed.
dr, A,
where
\j =U.001XJ (^^)
The forward difference method was seiected because it oniy required two calculations of
the cost function. This helped to speed up the time to convergence. The central
difference method wouid have required three calculations for each design variable. Ail of
the constraints are linear so their gradients were easily obtainable analytically. The final
calculation at this step is the maximum constraint violation V\(see eqn. 4.21).
The third step is to use the information from the first two steps to define the QP
subprobiem (eqns. 4. i 5-4. 17). At this point the QP subprobiem can be used to soive for
the search direction, d and Lagrange multipliers, u. In order to obtain these vaiues the
cost function, f(x) in the Lagrange Equation must be replaced by the QP cost function
(eqn. 4. 1 5) and the system of equations, 4. 1 1 -
4. 14, must be solved simultaneously.
The forth step is to check the convergence criteria to see that:
irn
-
*i
82
and
T7 _, _
If these conditions are satisfied the aigorithm has converged and the anaiysis can stop.
Otherwise continue to step five of the aigorithm.
kfe
The next step of the analysis is to modify the penalty parameter, R. For the
iteration the new penalty parameter Rk+i becomes:
* =
maxfo,rj (4.29)
where Rk is the existing penalty parameter and rk is the sum of the Lagrange multipliers
calculated in third step of the aigorithm. By updating the penaity parameter the necessary
condition wiii always be satisfied.
Next the step size must be determined. The inexact iine search method previously
developed was used here to calculate the proper step size. Once the step size is
determined the design point can be indexed. Tnerefore:
The finai step to the aigorithm is to index the counter, k=k+i, and repeat ail but
the first step. Iterations wiii continue untii convergence is reached.
4.3 Matiab Solution:
The CSD aigorithm was implemented for the optimization of the spindie shaft
using Matiab. The optimization applied to both the static and the dynamic models
developed in Chapters 3 and 4. The optimization program wiii return optimal vaiues for
the lateral stiffness and location of the spindie support bearings. As in previous chapters
83
the Matiab program reads in a batch fiie. The batch fiie aiiows the user to define the
spindie parameters. The optimization constants used by the CSD aigorithm were hard
coded into the program. Therefore the user of the program does not have the flexibility
to change these. The programming code used to perform the optimization anaiysis can be
found in Appendix D.
The sample spindie analyzed staticaiiy and dynamically in Chapters 2 and 3 was
aiso optimized to demonstrate the Optimization program. Figure 4.2 illustrates the batch
file read in by the program. The optimum parameters returned by the program for the
static and dynamic anaiysis are listed in tables 4.2 and 4.3 respectively.
Table 4.2 Optimum Values for Static Analysis
i Design Variable Optimum Value i Original Value

a_2 (m.) ! 7.47 j 7.5 I
a_3 (in.) 16.09 13.5
K_f(rbFm.) i 1000000 500000
K_r (Ibffin.) 1000000 100000
Table 4.3 Optimum Vaiues for Dynamic Analysis
Design Variable i Optimum Value ! Original Value i

! 7.3 7.5
j a_2(in.) ;
! a_3 Cm.) i 16.15 ! 13.5
! | 1000000 500000
K_f0b81n.)
K_r (Mln.) 1000000 100000
The static deflection of the existing spindie was approximateiy The
corresponding static deflection of the optimized spindie was approximateiy .00079". The
optimization reduced the deflection by a factor of 6. For the dynamic anaiysis the initial
deflection was about The corresponding deflection of the optimized spindle was
approximateiy .00025". Here the deflection improved by a factor of 9. The optimization
84
Batch File:
Geometry:
Number of Sections(#): _6J
Secton Leneth Outer Diameter Inner Diameter Area Moment ot Inertia

(#) (in) On) (in -21 fin-4-)
1 3 2.25 2 0.83449 0 47265783
2 3 2.375 2 125 0.58357 0.560861726
3 3 2.5 2.25 0.93265 0.659419991
4 3 2.625 2.375 0.98175 0.76890787
5 3 2 75 2.5 1.03084 0.889900605
6 3 2875 2 625 107992 1.022973438
~t r\
1 0
8 0 0
9 0 0
10 0 0
11 0 0
n
12 0
13 0 0
14 0 0
15 0 0
Bearings:
Lat Stiffness of RearBeanng (lb/In): 100000

Lat Stiffness of Front Bearing (lb/In ): 500000
Tor. Stiffness of Front Bearing (in- lb): 10000
Fraction of mom. on Front Bearing:
Locator! of Rear Bea-np** (i.n.) 7.5
Fuilev:
Location ot Pulley (in ~) 4.5

Mass of Pulley (lb-s"2/1n): 0 007763975
Static Belt Tension (lb): 50
h'armenie Drive Force (lb): 50
Drive Fnpnugnru f_H2\ 13333
Pulley Unbalance (Ib-s^): 0.002279727
Tool:
MaSS of Tool (ifrs'2/rn)'. v.v^jo/aai/
State Cutting Force (lb): 300

Harmonic Cutting Force (tb): 300
Cutting Frequency (Hz): 133.33
Tool Unbalance (Ib-s^): 0.004559453
Length of Tool (in)
Speed:
Spindle Shaft Speed (rpm): I 16.66666667
flrlatefiai Properties:

Density (lb/1n"3): 0.289
Figure 4.2 Batch Fiie for Sample Spindle Problem
85
improved the deflection dynamic
of the spindie significantly for both the static and
modeis.
86
5.0 Conclusion:
The anaiysis of a machine tooi spindie began by developing a modei to soive for
the static lateral deflection. A program was developed using Matlab that reads in the
geometry and ioad information for a spindie and reports piots of the iaterai deflection.
The user must simply enter the appropriate information into a spreadsheet, termed a
fiie"
"batch and the program wiii suppiy a piot of the spindies deformed shape. A sample
spindle was analyzed using this program and the FEA program Ansys. Both analyses
yielded comparable results.
Next the anaiysis was extended to include the dynamic response of the spindie.
Again a program was developed that wouid read in a batch fiie containing the appropriate
geometry and ioad information for the spindle. The program wouid then report piots of
the first four mode shapes of the spindie as weii as a piot of its dynamic forced response.
A sample spindie was aiso analyzed using the dynamic program and the results compared
to an FEA analysis. The FEA anaiysis agreed very closely with the dynamic analysis
program for both the mode shapes and iaterai deflection.
Finaiiy a program was developed that wouid optimize the spindie by minimizing
the deflection at the interface between the cutting tooi and the spindie shaft. The program
performs an optimization for both the static and dynamic analyses. The program
optimizes the location and stiffness of the spindie support bearings. A sample spindie
was optimized both statically and dynamically using the program. One key result of
these analyses was that the optimum parameters for both the static and dynamic anaiysis
were approximateiy the same. This is very important because the complexity and time
required to analyze and optimize the spindie than that
dynamically was significantly more
of the static anaiysis and optimization. Tnerefore if the goal of a spindie designer is to
optimize an existing spindie design, it couid be done statically with much iess effort and
time than it couid be done dynamically.
Although these programs are very powerful in designing a machine tooi spindie,
further work couid make them of more vaiue to a spindle designer. Tne first
recommendation wouid be to equip these programs with a graphical user interface. This
would make the interface between the designer and the program much more user friendly.
The next recommendation wouid be to make the programming code more
efficient. Currently it takes a considerable amount of time to process the dynamic
anaiysis and optimization. There are several iterations in each of these analyses. The
dynamic optimization takes several hours to run.
Another recommendation for future work wouid be to enhance the dynamic ioads
applied to the spindie modei. Currently harmonic ioads are assumed for the cutting force
and drive force. The accuracy of the forced response couid be further refined by taking
measurements of the cutting force and drive force for an existing machine tool spindles.
This couid be taken one step further by creating a database of the cutting forces for
various cutting toois. If a spindie was to primarily use one type of cutting tooi the ioad
information couid be read into the dynamic program from the database.
One final suggestion for future work would be to work with a

bearing
manufacturer to create a database of bearings and their stiffness vaiues. This wouid
88
provide the stiffness for the static and dynamic analyses as weii as the maximum
aiiowabie stiffness vaiues for the optimization analyses.
89
REFERENCES
[1] Ai-Shareef, K.J.H, Brandon, J. A., "On the Applicability ofModai and Response
Representations in the Dynamic Analysis ofMachine Tooi Spindie Bearing
Systems,"
Journal of Engineering Manufacture, Voi. 205, 1991, pp. 139-145.
[2] Ai-Shareef, K.J.H, Brandon, J. A., "On the Quasi-Static Design ofMachine Tooi
"
Spindies, Journal of Engineering Manufacture, Voi. 204, 1990, pp. 91-104.
[3] Arora, J, Introduction to Optimum Design, McGraw-Hiii Inc., Boston, 1989.
[4] Curti,G., Raffa, FA, Vatta, F., "Steady-State Unbalance Response of Continuous
"
Rotors on Anisotropic Supports, Vibration ofRotating Systems, Voi. 60, ASME

1993, pp. 27-34.
Prentice-
[5] Dahieh, M.D., Thomson, W.T., Theory of Vibration with Applications,
Hall, Inc., London, i99S.
"
[6] Lewinchai, L., "Machine Tooi Spindie Applications, SKF Industries, Inc.
Engineering and Research, SKF Norden, Feb 1983.
"
[7] Lewinchai, L., "Machine Tooi Spindie Applications, SKF Industries, Inc.
Engineering and Research, SKF Norden, Feb 1985.
[8] Montusiewicz,J ., Osycska,A, "Computer Aided Optimum Design ofMachine

Bearings"
Tool Spindie Systems with Hydrostatic Journal ofEngineering

Manufacture, Vol. 211, 1997, pp. 43-51.
[9] Murphy ,B., "Improved Rotordynamics Unbalance Response Calculations Using

the Polynomial Method", Vibration ofRotating Systems, Vol. 60, ASME 1993,
pp. 35-42.
90
Appendix A
Template for Batch Fiie
91
oaten r ne.
Geometry'.
Number of Sections(#): _6J

m (in) On) (in) (in"2) (irT4)
1 3 225 2 0.83449 0.47265783
2 3 2.375 2 125 0 88357 0 560861726
3 3 2.5 2.25 Q.S3255 0.659419991
4 3 2.625 2 375 0.88175 0 76890787
5 3 2 75 25 103084 0 889900605
6 3 2.875 2625 107992 1 022973438
7 0 0
8 0 0
9 0 0
n n
10
11 0 0
12 0 0
13 0 0
14 0 0
15 0 0
Bearings'
Let Stiffness of R ear Bearing (lb/in): 100000

Lat. Stiffness of Front Bearing (Ib^n.): 500000
Tor. Stiffness of Front Bearing (in-lb): 10000
Fraction of mom. on Front Bearing: 0.1
Location of Rear Bearing (in.)
rst Bearing {
Pulley:
Location of Pulley (in.) 4.5

Mass of Pulley (Ib-s^/ln): 0.007763975
Static Belt Tension (lb): 50
Harmonic Drive Force (lb): 50
Drive Frequency \Hzf. 133.33
D.rllp./ I InKalsnra /IK.S~?V
0002279727
Tool:
Mass of Tooi (ifa-s'2/in): 0.0258 /y91/

StaticCutting Force (lb): 300
Harmonic Cutting Force (lb): 300
Cutting Frequency (Hz): \3533
iGui unuaiafH.r \1u-3 .),
Length of TooJ (in)

Speed:
Spindle Shaft Speed (rpm'): 1 16.66666667 1
Modulus of Elasticity (psf): 30000000 ;

Density (lb/in~3): 0.289
92
Appendix B
Matiab Programs for the Static Anaiysis
93
% Filename : spindiec
% This fiie is the parent to ail other subroutines that wiii be employed
% to perform the static and dynamic analysis.

*C*
% Revision
dear
information from "Batch!.

wki"
% Reads in geometry and loading
constantsc
% Calculates the static deformation.
deformationc
% Calculates the dynamic response.
dynamicc
94
"Constantsc*
% Filename :
Vo This subroutine coilects the constants required to perform the static and dynamic
analysis.
"batch" "RevC"
% The subroutine requires a spreadsheet to be in the foider

"batch"
% contains aii the inputs required to drive the anaiysis.
dear;
temp
=
0;
=
ct wkIread(!Batehf);
=
metooi ct(49,4);
=
mepuiiey ct(4i,4);
=
speed ct(53,4);
tlunbai =
metool*(speed*2*pi)A2;
=
piunbai mepuliey*(speed*2*pi)A2;
ftharm =
ct(47,4);
fpharm =
ct(39,4);
freqpuiiey =
ct(40,4);
freqtooi =
ct(48,4);
E =
ct(57,4);
density =
ct(58,4);
N =
ct(5,4);
Mpuiley =ct(37,4);
Mtooi =
ct(45,4);
Kr =
ct(27,4);
Kf =
ct(28,4);
95
Ktf =
ct(29,4);
=
xl
ct(30,4);
tl =
ct(50,4);
fl=ct(38,4);
f2 =
ct(46,4);
fora=l:N;
b(n) =
ct(8+n,2)+temp;
temp
=
b(n);
Area(n) =
ct(8+n,5);
Inertia(n) =
ct(8+n,6);
R_out(n) =
ct(S+nJ)/2;
R_in(n) =
ct(8+n,4)/2;
end
Omsft =
ct(53,4);
=
a(l) ct(36,4);
=
a(2) ct(31,4);
=
a(3) et(32,4);
a(4)
=
b(N);
ifa(I)<=b(l)
XsecA(l) =
Area(l);
XsecI(l) =
Inertia(l);
XsecR(l) =
R_out(l);
96
end
for n=
1:N
ifb(n)<a(l)
ifb(n+l)>a(l)
XsecA(l) =
Area(n+1);
Xsecl( 1 ) =
Inertia(n+ 1 );
XsecR(l) =
R_out(n+l);
end
end
ifb(n)<a(2)
ifb(n+l)>a(2)
XsecA(2) =
Area(n+l);
Xsecl(2) =
Inertia(n+1);
XsecR(2) =
R_out(n+l);
end
end
ifb(n)<a(3)
ifb(n+l)>a(3)
XsecA(3) =
Area(n+1);
Xsecl(3) =
Inertia(n+1);
XsecR(3) =
R_out(n+l);
97
end
end
end
XsecA(4) =
Area(N);
Xsecl(4) =
Inertia(N);
XsecR(4) =
R_out(N);
forn=l:(N-l);
=
r(n) Inertia(n+1)/Inertia(n);
R(n) =
(Area(n+l)/Inertia(n+l))A.25/(Area(n)/Inertia(n))A.25;
end
x(l)
=
0;
forn=l:N
z =
2*n;
x(z)
=
b(n);
ifz<2*N;
x(z+l)
=
b(n);
end
end
forn=
1:N
z =
2*n-l,
Outer<z) =
R_out(n);
Outer(z+l) =
R_out(n);
98
Inner(z) =
R_in(n);
Inner(z+1)=
R_in(n);
end
plot(x,Outer,x,Inner,a(l),XsecR(l),,s,,a(2),XsecR(2),,*,,a(3),XsecR(3),,*',a(4),XsecR(4),,p
);
xlabel('x,(in)')
ylabel('shaft radius, (in)')
titleCinput dimensions')
legend('OD','ID','Pulley7Rear Bearing'/Front BearingVTool',2)
axis([0,b(N),0,(b(N))/2])
99
% Filename : deformationc
% This subroutine will calculate the static lateral deflection of the spindle
% Revision "C"
% Contribution due to bearing deformation
deflbrgc
% Contribution due to shaft deformation
deflshftc
yb +
=
yt ys;
% Creates Plots of the shaft deformation
plotsc
100
% Filename .
deflbrgc
% This subroutine will calcuate the deformation of the bearings
"C"
% Revision
x
=
linspace(0,(a(4)),100);
*
mb =
f2 ((a(4)+tl)-a(3))*xl;
Rl =
(fl*(a(3)-a(l))+mb-f2*((a(4)+tl)-a(3)))/(a(3)-a(2));
R2 =
fl +f2-Rl;
deltal =
-Rl/Kr;
delta2 =
-R2/Kf;
m
=
(delta2-deltal)/(a(3)-a(2));
for n=
1.100
yb(n)
=
m*(x(n)-a(2)) + deltal;
end
101
% Filename: deflshftc
% This subroutine will calculate the deformation of the elastic shaft
% in rigid supports
% This portion of the routine will scale the loads:
"C"
% Revision
for j= 1:3
n(j)
=
i;
for i =
2:N
ifa(j)>b(i -1)
ifa(j)<=
b(i)
=
n(j) ;
else
end
else
n(j)=n(j);
end
end
end
fle =
Inertia(N)/Inertia(n(l))*fl ;
Rle =
Inertia(N)/Inertia(n(2))*Rl;
mbe
=
Inertia(N)/Inertia(n(3))*mb;
R2e =
Inertia(N)/Inertia(n(3))*R2;
for i=l:(N-l)
102
Rl*sing(b(i),a(2),0) + R2*sing(b(i),a(3),0) fl*sing(b(i),a(l),0);
=
voi -
* *
moi =
Rl (b(i)-a(2)) sing(b(i),a(2),0) + R2*(b(i)-a(3))*sing(b(i),a<3),0).
-
fl*(Ki)-a(D)*sing(b(i),a(l),0) mb*sing(b(i),a(3),0); -
Inertia(N)*(l/Inertia(i)
=
v(i)
-
l/Inertia(i+l))*voi;
m(i)
=
Inertia(N)*(l/Inertia(i) -
l/Inertia(i+l))*moi;
end
x=
linspace(0,a(4),100);
tempi =0;
temp2 =
0;
forr=
1:(N-1)
tempi =
tempi + v(r)/6*sing(a(3),b(r),3) -
v(r)/6*sing(a(2),b(r),3)
+
m(r)/2*sing(a(3),b(r),2) m(r)/2*sing(a(2),b(r),2);
-
temp2 =
temp2 + v(r)/6*sing(a(2),b(r),3) + m(r)/2*sing(a(2),b(r),2);
end
=
+ tempi
ql l/(a(2)-a(3))*(fle/6*sing(a(3),a(l),3) Rle/6*sing(a(3),a(2),3) -
fle/6*sing(a(2),a(l),3));
=
q2 -fle/6*sing(a(2),a(l),3)
-
temp2 -
ql*a(2);
for r=
1:100
temp3 =
0;
for s=l:(N-l)
temp3 =
temp3 + v(s)/6*sing(x(r),b(s),3) + m(s)/2*sing(x(r),b(s),2);
end
ys(r)
=
-l/(E*Inertia(N))*(fle/6*sing(x(r),a(l),3)
-
Rle/6*sing(x(r),a(2),3).
103
-R2e/6*sing(x(r)X3),3) +temp3 +
mbe/2*sing(x(r),a(3),2) + ql*x(r) + q2);
end
104
% Filename plotsc
% This subroutine will create plots of the shaft deformation

"C"
% Revision
figure
plot(x,yb)
titleCDeflection Contribution of Bearings1)
figure
plot(x,ys)
title('Deflection Contribution of Elastic Shaft')
figure
plot(x,yt)
title('Combined Spindle Deflection')
105
function sing =
f(o,p,q)
% Function : Singularity
% This subroutine defines a new function to be used in subsequent calculations
if o<p
sing
=
0;
else
sing
=
(o-p)Aq;
end
106
Appendix C
Matlab Programs for the Dynamic Analysis
107
"dynamicc"
% Filename :
% Rev C
% This subroutine solves for the dynamic response of the spindle.
% Must read in the batch file prior to executing this subroutine.
% Solves for the eigenvalues.
frequencyc
% Solves for the eigenvectors.
modec
% Solves for the dynamic response.
forcedc
108
"frequencyc"
% Filename :
% This subroutine is used to plot the frequency equation.
% There needs to be a matrix A that drives this subroutine.
Omega =
linspace( 1,2000,50);
for n=l:50;
w
=
Omega(n);
XI =
(density*XsecA(l)*wA2/E/XsecI(l))A25*a(l);
X2 =
(density*XsecA(2)*wA2/EyXsed(2))A.25*a(2),
X3 =
(density*XsecA(3)*wA2/E/XsecI(3))A.25*a(3);
X4 =
Transferc
CI =
Mpulley*Xl/(density*XsecA(l)*a(l));
C2 =
Kr/(X2/a(2))A3/E/XsecI(2);
C3 =Ktf/(X3/a(3))A2/E/XsecI(3);
C4 =
Kf/(X3/a(3))A3/E/XsecI(3);
C5 =
Mtool*X4/(density*XsecA(4)*a(4));
templ(l,:)
=
[10-1 Oj;
templ(2,:)
=
[0 10-l];
Ql=templ*Tl;
temp2(l,:)
=
[-(sinh(Xl)-Cl*cosh(Xl)) -(cosh(Xl)-Cl*sinh(Xl))
-(sin(Xl)-
C 1 *cos(Xl )) (cos(Xl )+C 1 *sin(Xl ))];
temp2(2,:)
=
[-cosh(Xl) -sinh(Xl) -cos(Xl) -sin(Xl)];
temp2(3,:)
=
[-sinh(Xl) -cosh(Xl) sin(Xl) -cos(Xl)];
109
temp2(4,:)
=
[-cosh(Xl) -sinh(Xl) cos(Xl) sin(Xl)];
Q2 =
temp2*T2;
temp3(l,:)
=
[ sinh(X2) cosh(X2) sin(X2) -cos(X2)];
temp3(2,:)
=
[-cosh(X2) -sinh(X2) -cos(X2) -sin(X2)];
temp3(3,:)
=
[-sinh(X2) -cosh(X2) sin(X2) -cos(X2)],
temp3(4,:)
=
[-cosh(X2) -sinh(X2) cos(X2) sin(X2)];
Q3=temp3*T3;
temp4(l,:)
=
temp4(2,:)
=
temp4(3,:)
=
[-sinh(X3) -cosh(X3) sin(X3) -cos(X3)];
temp4(4,:)
=
[cosh(X3) sinh(X3) -cos(X3) -sin(X3)];
Q4 =
temp4*T4;
A(l,.) =
[Ql(l,l) Ql(l,2) Ql(l,3) Ql(l,4) 00000000000 0];
A(2,:) =
[Ql(2,l) Ql(2,2) Ql(2,3) Ql(2,4) 00000000000 0];
A(3,:) [sinh(Xl) cosh(Xl) Q2(l,l) Q2(l,2) Q2(l,3) Q2(l,4) 0 0 00

=
sin(Xl) -cos(Xl)
0 0 0 0];
A(4,:) =
[cosh(Xl) sinh(Xl) cos(Xl) sin(Xl) Q2(2,l) Q2(2,2) Q2(2,3) Q2(2,4) 0 00 0
0 0 0 0];
A(5,:) [sinh(Xl) Q2(3,l) Q2(3,2) Q2(3,3) Q2(3,4) 0 0 00

=
cosh(Xl) -sin(Xl) cos(Xl)
00 0 0];
A(6,:) [cosh(Xl) sinh(Xl)

=
-cos(Xl) -sin(Xl) Q2(4,l) Q2(4,2) Q2(4,3) Q2(4,4) 0 0 0
0 0 0 0 0];
A(7,:) (C2*cosh(X2)-sinh(X2)) (C2*sinh(X2)-cosh(X2)) (C2*cos(X2)-

[0000
=
sin(X2)) (C2*sin{X2)+cos(X2)) Q3(l,l) Q3(l,2) Q3(l,3) Q3(l,4) 0 0 0 0];
A(8,:) =
[0000 cosh(X2) sinh(X2) cos(X2) sin(X2) Q3(2,l) Q3(2,2) Q3(2,3) Q3(2,4)
110
0 0 00];
A(9,:) =
[0000 sinh(X2) cosh(X2) -sin(X2) cos(X2) Q3(3,l) Q3(3,2) Q3(3,3) Q3(3,4)
0 0 0 0];
A(10,:) [0000 cosh(X2) sinh(X2)

=
-cos(X2) -sin(X2) Q3(4,l) Q3(4,2) Q3(4,3)
Q3(4,4) 0 0 0 0];
A(l 1,;) =
[0 0000000
(C4*cosh(X3)-sinh(X3)) (C4*sinh(X3)-cosh(X3))
(C4*cos(X3)-sin(X3)) (C4*sin(X3)+cos(X3)) Q4(l,l) Q4(l,2) Q4(l,3) Q4(l,4)];
A(12,:) =
[0 0000000 cosh(X3) sinh(X3) cos(X3) sin(X3) Q4{2,1) Q4(2,2) Q4(2,3)
Q4(2,4)];
A(13,:) =
[0 0000000 sinh(X3) cosh(X3) -sin(X3) cos(X3) Q4(3,l) Q4(3,2) Q4(3,3)
Q4(3,4)];
A(14,:) [0 0 0 0 0 0 0 0 (C3*sinh(X3)-cosh(X3)) (C3*cosh(X3)-sinh(X3)) (-

=
C3*sin(X3)+cos(X3)) (C3*cos(X3)+sin(X3)) Q4(4,l) Q4(4,2) Q4(4,3) Q4{4,4)];
A{1 5,:) =
[0 00000000000 cosh(X4) sinh(X4) -cos(X4) -sin(X4)];
A(16,:) =
[0 00000000000 (sinh(X4)+C5*cosh(X4)) (cosh(X4)+C5*sinh(X4))
(sin(X4)+C5 *cos(X4)) (-cos(X4)+C5*sin(X4))];
fork =1:15
A(k,:)=A(k,:)/A(k,k);
for p =
k+l:16
A(p,:) =
A(p,:)-A(k,:)*A(p,k);
end
end
freq(n) =
A(16,16);
end
n=l;
for s =
2:50;
ill
freq(s)*freq(s-l);
=
sign
if sign < 0
if abs((freq(s)-freq(s-l ))/(Omega(s)-Omega(s- 1 )))<25
bisectc
end
end
end
112
"transferj"
% Filename :
% This subroutine creates the matrices to transform the constants.
Tl=eye(4);
T2=eye(4);
T3=eye(4);
T4=eye(4);
for counter =
1:(N-1)
betal =
(Area(counter)*density*wA2/(E*Inertia(counter)))A25,
beta2 =
(Area(counter+l)*density*wA2/(E*Inertia(counter+l)))A.25;
=
rl r(counter);
Rl =
R(counter);
11 =
b(counter);
flag =
0;
T =
zeros(4,4);
T(l,l) ((rl*RlA2+l)/2*(cosh(betal*ll)*cosh(beta2*ll)-
=
Rl*sinh(betal*ll)*sinh(beta2*ll)));
T(U) ((rl*RlA2+l)/2*(cosh(betal*ll)*sinh(beta2*ll)
= -
Rl *sinh(betal *ll)*cosh(beta2*ll)));
T(l,3) = -((rl*RlA2-l)/2*(cosh(betal*ll)*cos(beta2*ll)+
Rl *sinh(betal *ll)*sin(beta2*ll)));
T(l,4) =
-((rl*RlA2-l)/2*(cosh(betal*ll)*sin(beta2*ll)
-
Rl *sinh(betal *ll)*cos(beta2*ll)));
((rl*RlA2+l)/2*(Rl*cosh(betal*ll)*sinh(beta2*ll)-
T(2,l)
=
sinh(betal *ll)*cosh(beta2*ll)));
T(2,2)
=
((rl*RlA2+l)/2*(Rl*cosh(betal*ll)*cosh(beta2*ll) -
sinh(betal*ll)*sinh(beta2*ll)));
113
T(2,3) =
((rl *R1 A2-1)/2*(R1 *cosh(betal *11 )*sin(beta2*ll)+
sinh(betal*ll)*cos(beta2*ll)));
T(2,4) =
((rl*RlA2-l)/2*(sinh(betal*ll)*sin(beta2*ll>
Rl *cosh(betal *ll)*cos(beta2*ll)));
T(3,l) ((rl*RlA2-l)/2*(Rl*sin(betal*ll)*sinh(beta2*ll)
= -
cos(betal *ll)*cosh(beta2*ll)));
T(3,2) ((rl*RlA2-l)/2*(Rl*sin(betal*ll)*cosh(beta2*ll)
= -
cos(betal *11 )*sinh(beta2*ll )));
((rl*RlA2+l)/2*(cos(betal*ll)*cos(beta2*ll)+
T(3,3) =
Rl*sin(betal*ll)*sin(beta2*ll)));
T(3,4) ((rl*RlA2+l)/2*(cos(betal*ll)*sin(beta2*ll)
= -
Rl*sin(betal*ll)*cos(beta2*ll)));
T(4,l) = -((rl*RlA2-l)/2*(Rl*cos(betal*ll)*sinh(beta2*ll)+
sin(betal *1 1 )*cosh(beta2*l 1 )));
T(4,2) = -((rl*RlA2-l)/2*(Rl*cos(betal*ll)*cosh(beta2*ll)+
sin(betal *ll)*sinh(beta2*ll)));
T(4^) =
((rl *R1 A2+l)/2*(sin(betal *ll)*cos(beta2*ll)-
Rl *cos(betal *ll)*sin(beta2*ll)));
((rl*RlA2+l)/2*(sin(betal*ll)*sin(beta2*ll)+
T(4,4) =
Rl *cos(betal *ll)*cos(beta2*ll)));
if a(3)<
b(counter)
T4 =
T4*T;
flag =
l;
end
if a(2)<
b(counter)
if flag <1;
T3 =
T3*T;
114
flag=l;
end
end
if a(l)<b(counter)
if flag <1;
T2 =
T2*T;
flag=l;
end
end
if flag < 1;
T1 =
T1*T;
flag=0;
end
end
115
"bisectc"
% Filename :
% This subroutine refines the incremental root finding search
% using the bisection method.
Xlow =
Omega(s-l);
Xhigh =
Omega(s);
freq_ref =
freq(s-l);
ffeq_new =
freq(s-l);
while abs(Xhigh-Xlow)>le-12;
Xnew =
.5*(Xlow+Xhigh);
w
=
Xnew;
XI =
(density*XsecA(l)*wA2/E/XsecI(l))A.25*a(l);
X2 =
(density*XsecA(2)*wA2/E/Xsed(2))A.25*a(2);
X3 =
X4 =
Transferc
CI =
C2 =Kr/(X2/a(2))A3/E/XsecI(2);
C3 =
Ktf/(X3/a(3))A2/E/XsecI(3);
C4 =
C5 =
templ(l,:)
=
[10-10];
templ(2,:)
=
[0 10-l];
Ql= tempi *T1;
116
temp2(l,:)
=
[-(sinh(Xl)-Cl*cosh(Xl)) -(cosh(Xl)-Cl*sinh(Xl)) -(sin(Xl)-
Cl*cos(Xl)) (cos(Xl)+Cl*sin(Xl))];
temp2(2,:)
=
temp2(3,:)
=
temp2(4,:)
=
Q2 =
temp2*T2;
temp3(l,:)
=
temp3(2,:)
=
temp3(3,:)
=
temp3(4,:)
=
f-cosh(X2) -sinh(X2) cos(X2) sin(X2)];
Q3=temp3*T3;
temp4(l,:)
=
temp4(2,:)
=
temp4(3,:)
=
temp4(4,:)
=
Q4 =
temp4*T4;
A(l,:) =
[Ql(l,l) Ql(l,2) Ql(l,3) Ql(l,4) 00000000000 0];
A(2,:) =
[Ql(2,l) Ql(2,2) Ql(2,3) Ql(2,4) 00000000000 0];
A(3,:) =
[sinh(Xl) cosh(Xl) sin(Xl) -cos(Xl) Q2(l,l) Q2(l,2) Q2(l,3) Q2(l,4) 0 000
0 0 0 0];
A(4,:) =
[cosh(Xl) sinh(Xl) cos(Xl) sin(Xl) Q2(2,l) Q2(2,2) Q2(2,3) Q2(2,4) 0 0 0 0
0 0 0 0];
A(5,:) =
[sinh(Xl) cosh(Xl) -sin(Xl) cos(Xl) Q2(3,l) Q2(3,2) Q2(3,3) Q2(3,4) 0000
0 0 0 0];
117
A(6,:)
=
[cosh(Xl) sinh(Xl) -cos(Xl) -sin(Xl) Q2(4,l) Q2(4,2) Q2(4,3) Q2(4,4) 0 0 0
0 0 00 0];
A(7,:) =
(C2*cosh(X2)-sinh(X2)) (C2*sinh(X2)-cosh(X2)) (C2*cos(X2)-
[0000
sin(X2)) (C2*sin(X2)+cos(X2)) Q3(l,l) Q3(l,2) Q3(l,3) Q3(l,4) 0 0 0 0];
A(8,:) =
0 0 0 0];
A(9,:) =
000 0];
A(10,:) [0000 cosh(X2)

=
sinh(X2) -cos(X2) -sin(X2) Q3(4,l) Q3(4,2) Q3(4,3)
Q3(4,4) 0 0 0 0];
A(l 1,:) =
[0 0000000
(C4*cos(X3)-sin(X3)) (C4*sin(X3)+cos(X3)) 04(1,1) Q4(l,2) Q4(l,3) 04(1,4)];
A(12,:) =
[0 0000000 cosh(X3) sinh(X3) cos(X3) sin(X3) Q4(2,l) Q4(2,2) Q4(2,3)
04(2,4)];
A(13,:) =
Q4(3,4)];
A(14,:) [0 0000000 (C3*sinh(X3)-cosh(X3)) (C3*cosh(X3)-sinh(X3)) (-

=
C3*sin(X3)+cos(X3)) (C3*cos(X3)+sin(X3)) Q4(4,l) Q4(4,2) Q4(4,3) Q4(4,4)];
A(15,:) =[0 00000000000 cosh(X4) sinh(X4) -cos(X4) -sin(X4)];
A(16,:) =
(sin(X4)+C5 *cos(X4)) (-cos(X4)+C5 *sin(X4))];
fork =1:15
A(k,:)=A(k,:)/A(k,k);
for p =
k+l:16
A(p,:) =
A(p,:)-A(k,:)*A(p,k);
end
end
118
freq_new =
A( 1 6, 1 6);
if (freq_new*freq_ref)< 0
Xhigh =
Xnew;
else
Xlow =
Xnew;
freq_ref =
freq_new;
end
end
Root(n) =
Xnew; % This will be in rad/s
Frequency(n) =
Xnew/2/pi;
=
n n+l;
119
"modec"
% Filename :
% This subroutine calculates the first (4) mode shapes.
globai c
for m =
1 :4
w =
Root(m);
XI =
(density*XsecA(l)*wA2/E/Xsed(l))A.25*a(l);
X2 =
X3 =
X4 =
(density *XsecA(4)*wA2/E/XsecI(4))A.25*a(4);
Transferc
CI =
C2 =
C3 =
Ktf/(X3/a(3))A2/E/XsecI(3);
C4 =
C5 =
templ(l,:)
=
[10-10];
templ(2,:)
=
[010 -1];
Ql =
templ*Tl;
temp2(l,:) [-(sinh(Xl)-Cl*cosh(Xl))
=
-(cosh(Xl)-Cl*sinh(Xl)) -(sin(Xl>
Cl*cos(Xl)) (cos(Xl)+Cl*sin(Xl))];
temp2(2,:)
=
temp2(3,:)
=
[-sinh(Xl) -cosh(Xl) sin(Xl) -cos(Xl)],
temp2(4,:)
=
120
Q2 =
temp2*T2;
temP3(l,:)
=
temp3(2,:)
=
temp3(3,:) =
temp3(4,:)
=
Q3=temp3*T3;
temp4(l,:)
=
temp4(2,:)
=
temp4(3,:)
=
temp4(4,:)
=
Q4=temp4*T4;
A(l,:) =
[Ql(l,l) Ql(l,2) Ql(l,3) Ql(l,4) 00000000000 0];
A(2,:) =
[Ql(2,l) Ql(2,2) Ql(2,3) Ql(2,4) 00000000000 0];
A(3,:) =
[sinh(Xl) cosh(Xl) sin(Xl) -cos(Xl) Q2(l,l) Q2(l,2) Q2(l,3) Q2(l,4) 00 0 0
0 0 0 0];
A(4,:) [cosh(Xl) Q2(2,l) Q2(2,2) Q2(2,3) Q2(2,4) 0 0 00

=
sinh(Xl) cos(Xl) sin(Xl)
0 0 0 0];
A(5,:) [sinh(Xl) cosh(Xl) Q2(3,l) Q2(3,2) Q2(3,3) Q2(3,4) 0 0 0 0

=
-sin(Xl) cos(Xl)
0 00 0];
A(6,:) [cosh(Xl) sinh(Xl)

=

0 0 0 0 0];
A(7,:) =
(C2*cosh(X2)-sinh(X2)) (C2*sinh(X2)-cosh(X2)) (C2*cos(X2>-
[0000
A(8,:) =
0 0 0 0];
121
A(9,:) =
0 00 0];
A(10,:) [0000 cosh(X2) sinh(X2)

=
-cos(X2) -sin(X2) Q3(4,l) Q3(4^) Q3(4,3)
Q3(4,4) 0 0 0 0];
A(l 1,:) =
[0 0 0 0 0 0 0 0
A(12,:) =
Q4(2,4)];
A(13,:) =
Q4(3,4)];

=
A(15,:) =
A(16,.) [0 00000000000 (sinh(X4)+C5*cosh(X4)) (cosh(X4)+C5*sinh(X4))

=
(sin(X4)+C5*cos(X4))(-cos(X4)+C5*sin(X4))];
fork =1:15
A(k,:)=A(k,:)/A(k,k);
forp=k+l:16
A(p,:) =
A(p,:)-A(k,:)*A(p,k);
end
end
c(16^n)=l;
fork =1:15;
g=16-k;
temp
=
0;
for p =
g+l:16;
122
= + A(g,p)*c(p,m);
temp temp
end
=
c(g,m) -temp;
end
xl =linspace(0,a(4),31);
for n=
1:31
=
x xl(n);
delta =
12;
if xl(n)<a(3)
delta =
8;
if xl(n)<a(2)
delta =
4;
ifxl(n)<a(l),
delta =
0;
end
end
end
z=l;
fori=l:N
ifx>=b(i);
z =
i;
end
123
end
=
ydyn(n,m)
c((delta+l),m)*cosh(x*(density*Area(z)*wA2/E/Inertia(z))A.25)+c((delta+2)^n)*sinh(x*(
density*
Area(z)*wA2/E/Inertia(z))A.25)+c((delta+3),m)*cos(x*(density*Area(z)*wA2/E/I
neitia(z))A.25)+c((deltii+4),m)*sin(x*(density*Area(z)*wA2/E/Inertia(z))A.25);
end
figure
plouxl^dynCm),'*')
L =
sprintf('mode %0.5g',m);
title(L)
end
124
"modec"
% Filename :
% This subroutine calculates the first (4) mode shapes.
global c
for m =
1 :4
w =
Root(m);
XI =
(density*XsecA(l)*wA2/E/XsecI(l))A25*a(l);
X2 =
X3 =
(density*XsecA(3)*wA2/E/XsecI(3))A.25*a(3),
X4 =
Transferc
CI =
C2 =
C3 =
Ktf7(X3/a(3))A2/E/XsecI(3);
C4 =
Kf7(X3/a(3))A3/E/XsecI(3);
C5 =
templ(l,:)
=
[10-10];
templ(2,:)
=
[0 10-l];
Ql= tempi *T1;
[-(sinh(Xl)-Cl*cosh(Xl)) -(cosh(Xl)-Cl*sinh(Xl)) -(sin(Xl>

=
temp2(l,:)
Cl*cos(Xl)) (cos(Xl)-rCl*sin(Xl))];
temp2(2,:)
=
temp2(3,:)
=
temp2(4,:)
=
125
Q2 =
temp2*T2;
temp3(l,:)
=
temp3(2,:)
=
temp3(3,:) =
temp3(4,:)
=
Q3=temp3*T3;
temp4(l,:)
=
temp4(2,:)
=
temp4(3,:)
=
temp4(4,:)
=
Q4 =
temp4*T4;
A(l,:) =
[Ql(l,l) Ql(l,2) Ql(l,3) Ql(l,4) 00000000000 0];
A(2,:) =
[Ql(2,l) Ql(2,2) Ql(2,3) Ql(2,4) 00000000000 0];
A(3,:) =
[sinh(Xl) cosh(Xl) sin(Xl) -cos(Xl) Q2(l,l) Q2(l,2) Q2(l,3) Q2(l,4) 0000
0 0 0 0];
A(4,:) [cosh(Xl) Q2(2,l) Q2(2,2) Q2(2,3) Q2(2,4) 0 0 00

=
sinh(Xl) cos(Xl) sin(Xl)
0 0 0 0];
A(5,:) [sinh(Xl) cosh(Xl) -sin(Xl) Q2(3,l) Q2(3,2) Q2(3,3) Q2(3,4) 0 0 0 0

=
cos(Xl)
000 0];
A(6,;) [cosh(Xl) sinh(Xl)

=
0 0 0 0 0];
A(7,:) =
(C2*cosh(X2)-sinh(X2)) (C2*sinh(X2)-cosh(X2)) (C2*cos(X2)-
[0000
A(8,:) =
0 0 0 0];
126
A(9,:) =
00 00];
A(10,:) [0000 cosh(X2)

=
sinh(X2) -cos(X2) -sin(X2) Q3(4,l) Q3(4,2) Q3(4,3)
Q3(4,4) 0 0 0 0];
A(l 1,:) =
[0 0 0 0 0 0 0 0
A(12,:) =
Q4(2,4)];
A(13,:) =
Q4(3,4)];

=
A(15,:) =
A(l 6,:) =
(sin(X4)+C5 *cos(X4)) (-cos(X4)+C5 *sin(X4))];
for k= 1:15
A(k,:)=A(k,:)/A(k,k);
for p =
k+1.16
A(p,:) =
A(p,:)-A(k,:)*A(p,k);
end
end
c(16,m)=l;
fork =1:15;
g=16-k;
temp
=
0;
for p =
g+l:16;
127
= + A(g,p)*c(p,m);
temp temp
end
=
c(g,m) -temp;
end
x1 =
linspace(0,a(4),3 1 );
for n=
1:31
=
x xl(n);
delta =
12;
if xl(n)<a(3)
delta =
8;
ifxl(n)<a(2)
delta =
4;
if xl(n)<a(l);
delta =
0;
end
end
end
z=l;
for i =
1 :N
ifx>=b(i);
z
=
i;
end
128
end
=
ydyn(n,m)
c((delta+l),m)*cosh(x*(density*Area(z)*wA2/E/Inertia(z))A.25)+c((delta+2),m)*sinh(x*(
density*
Area(z)*wA2/E/Inertia(z))A.25)+c((delta+3),m)*cos(x*(density*Area(z)*wA2/E/I
nertia(z))A.25)+c((delta+4),m)*sin(x*(density*Area(z)*wA2/E/Inertia(z))A.25);
end
figure
plot(xl,ydyn(:,m),'*')
L =
sprintf^'mode %0.5g',m);
title(L)
end
129
function integ =
f(XH,XL,md,sec,beta)
% Function : integral of phiA2 from Xlow to Xhigh
% This subroutine defines a new function to be used in subsequent calculations
global c
integ =
lête*(c<sec,md)A2/2*(cosh(beta*XH)*sinh(beta*XH)+oeta*XH) ...
+ c(sec+l ,md)A2/2*(cosh(beta*XH)*sinh(beta*XH)-beta*XH).
. .
+
c(sec+2,md)A2/2*(cos(beta*XH)*sin(beta*XH)+beta*XH) . .
+ c(sec+3,md)A2/2*(-cos(beta*XH)*sin(beta*XH)+beta*XH). .
+ c(sec,md)*c(sec+l,md)*(sinh(beta*XH))A2 +
c(sec,md)*c<sec+2,nid)*(cosh(bete*XH)*sin(beta*XH)înh(beta*XH)*cos(beta*XH)) .
+ c(sec,md)*c(sec+3,md)*(sinh(beta*XH)*sin(beta*XH)-
cosh(beta*XH)*cos(beta*XH)). .
c<sec+l,nid)*c(sec+2,md)*(sinh(beta*XH)*sin(beta*XH)+cosh(beta*XH)*cos(beta*XH)
)...
+ c(sec+l,md)*c(sec+3,md)*(cosh(beta*XH)*sin(beta*XH)-
sinh(beta*XH)*cos(beta*XH))..
+
c(sec+2,md)*c(sec+3,md)*(sin(beta*XH))A2)
-
l/êta*(c(sec,md)A2/2*(cosh(beta*XL)*sinh(beta*XL)+oeta*XL)
+ c(sec+l,md)A2/2*(cosh(beta*XL)*sinh(beta*XL)-beta*XL).
+ c(sec+2,md)A2/2*(cos(beta*XL)*sin(beta*XL)+beta*XL). .
+ c(sec+3,md)A2/2*(-cos(beta*XL)*sin(beta*XL)4beta*XL). . .
+ c(sec,md)*c(sec+l,md)*(sinh(beta*XL))A2 +
c(sec,md)*c(sec+2,md)*(cosh(beta*XL)*sin(beta*XL)+sinh(beta*XL)*cos(beta*XL))..
+ c(sec,md)*c(sec+3,md)*(sinh(beta*XL)*sin(beta*XL)-
cosh(beta*XL)*cos(beta*XL)). .
o(sec+l,md)*c(sec+2,md)*(sinh(beta*XL)*sin(beta*XL)+cosh(beta*XL)*cos(beta*XL))
+ c(sec+l4nd)*c(sec+3^nd)*(cosh(beta*XL)*sin(beta*XL)-
sinh(beta*XL)*cos(beta*XL)). . .
+ c(sec+2,md)*c(sec+3,md)*(sin(beta*XL))A2);
130
Appendix D
Matlab Programs for the Optimization
131
% filename : "optimize"
% this file launches the optimization subroutines
constantso
stop
=
0;
for nn=
1:1000;
if stop <=
0;
%Step 1;
kk=l;
R(l)=l;
=
gamma .5;
epsl =.01;
=
eps2 .01;
Hess =
eye(4);
Krmax =
le6;
Kf_max=
le6;
DELTA =5;
OH =
2;
=
xn(l) a(2)/(a(l)+DELTA);
=
xn(2) a(3)/(a(4)-OH);
xn(3)
=
Kf/Kfmax;
xn(4)
=
Kr/Kr_max;
% Step 2
132
deformationo
yt(100)
constraint
=
temp max(G);
if temp <0
V =
0;
else
V =
temp;
end
gradient
dytn(l) =
dyt(l)*(a(l)+DELTA)
dytn(2) =
dyt(2)*(a(4)-OH)
dytn(3) =
dyt(3)*Kf_max
dytn(4) =
dyt(4)*Kr_max
% Step 3
multiplier
% Step 4
if abs_dd<epsl
stop=
1;
if V > eps2
stop
=
0;
133
end
end
% Step 5
+
uu(2) + uu(3)
=
rr
uu(l) +
uu(4);
R((kk+1)) =
max(rT,R(kk));
% Step 6
step_size
% Step 7
kk =
kk+l;
end
end
134
"gradient"
% filename
% this subroutine calculates the gradient of yt with respect to the design variables
% the design variables are Kf, Kr, a(l), a(2), a(3)
deformationo
=
yt_k yt(100);
% the gradient of yt w.r.t. a(2):
Deltal =
a(2)*001;
=
a(2) a(2)+Deltal;
deformationo
yt_kp yt(100);
dyt(l) =
(yt_kp-yt_k)/Deltal;
% the gradient of yt w.r.t. a(3):
Delta2 =
a(3)*01;
=
a(3) a(3)+Delta2;
deformationo
=
yt_kp yt(100);
dyt(2) =
(yt_kp-yt_k)/Delta2;
% the gradient of yt w.r.t. Kf:
Delta3=Kf*01;
Kf=Kf+Delta3;
deformationo
=
yt_kp yt(100);
dyt(3) =
135
% the gradient of yt w.r.t. Kr:
Delta4=Kr*.01;
Kr =
Kr+Delta4;
deformationo
=
yt_kp yt(100);
dyt(4) =
136
"constraint"
% Filename :
% This subroutine calculates the constraint violations and their derivatives.
G =
[(-xn(l)+l) (xn(2>l) (xn(3)-l) (xn(4)-l)];
dGl =
[-10 0 0];
dG2 =
[0 10 0];
dG3 =
[001 0];
dG4 =
[0 0 0 1];
137
% Filename "multiplier" .
% This file calculates the Lagrange Multipliers
clear AACC
if xn(l)>=l
uu(l)
=
0;
=
ss(l) sqrt(xn(l)-l);
swl =
0;
else
ss(l)
=
0;
=
swl 1;
end
if xn(2) <=
1
uu(2)
=
0;
=
ss(2) sqrt(xn(2)-l);
sw2 =
0;
else
ss(2)
=
0;
sw2=
1;
end
if xn(3) <=
1
uu(3)
=
0;
=
sw3 =
0;
138
else
ss(3)
=
0;
=
sw3 1;
end
if xn(4) <=
1
uu(4)
=
0;
sw4 =
0;
=
else
ss(4)
=
0;
sw4 =
1;
end
binary =
swl *2A3+sw2*2A2+sw3*2Al+sw4*2A0
if binary =
0
multO
end
if binary =
1
multl
end
if binary =
2
mult2
end
139
if binary =
2
mult2
end
if binary =
3
mult3
end
if binary ==
4
mult4
end
if binary =
5
mult5
end
if binary =
6
muh6
end
if binary =
7
mult7
end
if binary =
8
mult8
end
if binary =
9
140
mult9
end
if binary ==
10
mult 10
end
if binary =
11
mult 11
end
if binary =
12
mult 12
end
if binary =
13
mult 13
end
if binary =
14
mult 14
end
if binary =
15
mult 15
end
=
abs_dd
sqrt(dd(l)A2+dd(2)A2+dd(3)A2+dd(4)A2)
%ifswl=0
141
% ss(l) =
sqrt(dd(l)-dd(2)-DELTA/a(4));
%end
%ifsw2 =
0
% ss(2) =
sqrt(dd(3)-l);
%end
%ifsw3 =
0
% ss(3) =
sqrt(dd(4)-l);
%end
%ifsw4 =
0
% ss(4) =
sqrt(dd(5)-l);
%end
142
"multiplier"
% Filename .
% This file calculates the Lagrange Multipliers
clear AA CC
ifxn(l)>=l
uu(l)
=
0;
=
ss(l) sqrt(xn(l)-l);
swl=0;
else
ss(l)
=
0;
swl =
1;
end
if xn(2) <=
1
uu(2)
=
0;
=
sw2 =
0;
else
ss(2)
=
0;
sw2 =
1;
end
if xn(3) <=
1
uu(3)
=
0;
=
sw3 =
0;
143
else
=
ss(3) 0;
=
sw3 1;
end
ifxn(4)<=l
uu(4)
=
0;
sw4 =
0;
=
else
ss(4)
=
0;
sw4=
1;
end
binary =
swl*2A3+sw2*2A2+sw3*2Al+sw4*2A0
if binary =
0
multO
end
if binary =
1
multl
end
if binary =
2
mult2
end
144
if binary =
2
mult2
end
if binary =
3
mult3
end
if binary ==
4
mult4
end
if binary =
5
mult5
end
if binary =
6
muh6
end
if binary =
7
mult7
end
if binary =
8
mult8
end
if binary =
9
145
mult9
end
if binary =
10
mult 10
end
if binary =11
mult 11
end
if binary =
12
mult 12
end
if binary =13
mult 13
end
if binary =
14
mult 14
end
if binary =
15
mult 15
end
=
abs_dd
sqrt(dd(l)A2+dd(2)A2+dd(3)A2+dd(4)A2)
%ifswl=0
146
% ss(l) =
sqrt(dd(l)-dd(2)-DELTA/a(4));
%end
%ifsw2 =
0
% ss(2) =
sqrt(dd(3)-l);
%end
%ifsw3 =
0
% ss(3) =
sqrt(dd(4)-l);
%end
%ifsw4 =
0
% ss(4) =
sqrt(dd(5)-l);
%end
147
"multO"
% Filname :
AA(1,:) [(Hess(l,l)) (Hess(l,2)+Hess(2,l)) (Hess(l,3)+Hess(3,l))

=
(Hess(l,4)+Hess(4,l))];
AA(2,:) [(Hess(2,l)+Hess(l,2)) (Hess(2,2)) (Hess(2,3)+Hess(3,2))

=
(Hess(2,4)+Hess(4,2))];
AA(3,:) [(Hess(3,l)+Hess(l,3)) (Hess(3,2)+Hess(2,3)) (Hess(3,3))

=
(Hess(3,4)+Hess(4,3))];
AA(4,:) [(Hess(4,l)+Hess(4,l)) (Hess(4,2)+Hess(2,4)) (Hess(4,3)+Hess(3,4))

=
(Hess(4,4))];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
dd =
inv(AA)*CC
148
% Filname :
"multO"
AA(1,:) =
[(Hess(l,l)) (Hess(l,2)+Hess(2,l)) (Hess(l,3)+Hess(3,l))

(Hess(l,4)+Hess(4,l))];
AA(2,:) =
[(Hess(2,l)+Hess(l,2)) (Hess(2,2)) (Hess(2,3)+Hess(3,2))

(Hess(2,4)+Hess(4,2))];
AA(3,:) =
[(Hess(3,l)+Hess(l,3)) (Hess(3,2)+Hess(2,3)) (Hess(3,3))

(Hess(3,4)+Hess(4,3))];

=
(Hess(4,4))];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
dd =
inv( AA)*CC
149
"mult2"
% Filname :
dd(4) =
-G(3);
AA(1,:) =
[(Hess(l,l)) (Hess(l,2)+Hess(2,l)) (Hess(l,4)+Hess(4,l)) 0];
AA(2,:) =
[(Hess(2,l)+Hess(l,2)) (Hess<2,2)) (Hess(2,4)+Hess(4,2)) 0];
AA(3,:) =
[(Hess(3,l)+Hess(l,3)) (Hess(3,2)+Hess(2,3)) (Hess(3,4)+Hess(4,3)) 1];
AA(4,:) =
[(Hess(4,l)+Hess(4,l)) (Hess(4,2)+Hess(2,4)) Hess(4,4) 0];
CC(1,.) =
[-dytn(l)-(Hess(l,3)+Hess(3,l))*dd(3)];
CC(2,:) =
[-dytn(2)-(Hess(2,3)+Hess(3,2))*dd(3)];
CC<3,:) =
[-dytn(3)-Hess(3,3)*dd(3)];
CC(4,:) =
[-dytn(4)-(Hess(3,4)+Hess(4,3))*dd(3)];
tempmult =
inv(AA)*CC;
dd( 1 1 ) temp_muh( 1 1 );
=
, ,
dd(2, 1 ) temp_mult(2, 1 );
=
dd(4,l) =
temp_mult(3,l);
u(3)
=
temp_mult(4, 1);
150
"mult3"
%Filname:

=
(Hess(l,4)+Hess(4,l)) 0 0];

=
(Hess(2,4)+Hess(4,2)) 0 0];

=
(Hess(3,4)+Hess(4,3)) 1 0];
AA(4,:) =
[(Hess(4,l)+Hess(4,l)) (Hess(4,2)+Hess(2,4)) (Hess(4,3)+Hess(3,4))
(Hess(4,4))0 1];
AA(5,:) =
[001000];
AA(6,:) =
[000100];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,.) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
CC(5,:) =
[-G(3)];
CC(6,:) =
[-G(4)];
temp_mult =
inv(AA)*CC
=
dd(l,l) temp_mult(l,l)
=
dd(2,l) temp_mult(2,l)
=
dd(3,l) temp_mult(3,l)
dd(4, 1 ) =
temp_mult(4, 1 )
=
u(3) temp_mult(5,l);
=
u(4) temp_mull(6,l);
151
"mult4"
% Filname :
dd(2,l) =
-G(2);
AA(1,:) =
[(Hess(l,l)) (Hess(l,3)+Hess(3,l)) (Hess(l,4)+Hess(4,l)) 0];
AA(2,:) =
[(Hess(3,l)+Hess(l,3)) (Hess(3,2)+Hess(2,3)) (Hess(3,4)+Hess(4,3)) 1];
AA(3,:)
=
[(Hess(3,l)+Hess(l,3)) (Hess(3,3)) (Hess(4,3)+Hess(3,4)) 0];
AA(4,:) =
[(Hess(4,l)+Hess(l,4)) (Hess(4,3)+Hess(3,4)) (Hess(4,4)) 0];
CC(1,:) =
[-dytn(l)-(Hess(l,2)+Hess(2,l))*dd(2)];
CC(2,:) =
[-dytn(2)-(Hess(2,2))*dd(2)];
CC(3,:) =
[-dytn(3HHess(3,2)+Hess(2,3))*dd(2)];
CC(4,:) =
[-dytn(4)-(Hess(4,2)+Hess(2,4))*dd(2)];
tempmult =
inv(AA)*CC;
dd( 1 1 ) temp_muh( 1 1 );
=
, ,
dd(3 1 ) temp_mult(2, 1 );
=
,
dd(4,l) =
temp_mult(3,l);
u(2)
=
temp_mult(4, 1 );
152
"multi"
% Filname :
AA(1,.) [(Hess(l,l)) (Hess(l,2)+Hess(2,l)) (Hess(l,3)+Hess(3,l))

=
(Hess(l,4)+Hess(4,l)) 0 0];

=
(Hess(2,4)+Hess(4,2)) 1 0];

=
(Hess(3,4)+Hess(4,3)) 0 0];
AA(4,:) =
(Hess(4,4))0 1];
AA(5,:) =
[010000];
AA(6,:) =
[000100];
CC(1,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
CC(5,:) =
[-G(2)];
CC(6,:) =
[-G(4)];
tempmult =
inv(AA)*CC;
dd( 1 1 ) temp_mult( 1 1 );
=
, ,
=
dd(2,l) temp_mult(2,l);
=
dd(4,l) =
temp_mult(4,l);
u(2)
=
temp_mult(5, 1);
=
153
"mult6"
% Filname :

=
(Hess(l,4)+Hess(4,l)) 0 0];

=
(Hess(2,4)+Hess(4,2)) 1 0];

=
(Hess(3,4)+Hess(4,3)) 0 1];

=
(Hess(4,4)) 0 0];
AA(5,:) =
[010000];
AA(6,:) =
[001000];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
CC(5,:) =
[-G(2)];
CC(6,:) =
[-G(3)];
tempmult =
inv(AA)*CC;
dd( 1 1 ) temp_mult( 1 1 );
=
, ,
dd(2,l) =
temp_mult(2,l);
dd(3 1 )
,
=
temp_mult(3 , 1 );
dd(4,l) =
temp_mult(4,l);
u(2) =temp_mult(5,l);
154
"mult7"
% Filname

=
(Hess(l,4)+Hess(4,l)) 0 0 0],

=
(Hess(2,4)+Hess(4,2)) 1 0 0];

=
(Hess(3,4)+Hess(4,3)) 0 1 0];

=
(Hess(4,4)) 0 0 1];
AA(5,:) =
[0100000];
AA(6,:) =
[0010000];
AA(7,:) =
[0001000];
CC(l,:) =
[-dytn(l)];
CC(2,.) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
CC(5,:) =
[-G(2)];
CC(6,:) =
[-G(3)];
CC(7,:)
=
[-G(4)];
tempmult =
inv(AA)*CC;
dd(l,l)
=
temp_mult(l,l);
dd(2,l) =
temp_mult(2,l);
=
dd(4,l) =
temp_mult(4,l);
155
=
u(2) temp_mult(6, 1 );
=
u(4) temp_mult(8, 1);
156
"mult8"
% Filname :

=
(Hess(l,4)+Hess(4,l)) 1];

=
(Hess(2,4)+Hess(4,2)) 0];

=
(Hess(3,4)+Hess(4,3)) 0];

=
(Hess(4,4))0];
AA(5,:) =
[10000];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
CC(6,:)
=
[-G(1)];
tempmult =
inv(AA)*CC;
dd( 1 1 ),
=
temp_mult( 1,1);
dd(2,l) =
temp_mult(2,l);
dd(3,l) =
temp_mult(3,l);
dd(4,l) =
temp_mult(4,l);
u( 1 ) temp_mult(5, 1 );
=
157
"mult9"
% Filname :

=
(Hess(l,4)+Hess(4,l)) 0]; -1

=
(Hess(2,4)+Hess(4,2)) 0 0];

=
(Hess(3,4)+Hess(4,3)) 0 0];
AA(4,:) =
(Hess(4,4))0 1];
AA(6,:) =
[-10000 0];
AA(7,:)
=
[000100];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
CC(5,:) =
[-G(1)];
CC(6,:) =
[-G(4)];
temp_mult =
inv(AA)*CC;
dd(l,l)=temp_mult(l,l);
dd(2,l) =
temp_mult(2,l);
dd(3,l) =
temp_mult(3,l);
dd(4,l) =
temp_mult(4,l);
u( 1 )
=
temp_mult(5, 1 );
=
158
10"
%Filname: "mult

=
(Hess(l,4)+Hess(4,l)) 0]; -1

=
(Hess(2,4)+Hess(4,2)) 0 0];

=
(Hess(3,4)+Hess(4,3)) 0 1];
AA(4,:) [(Hess<4,l)+Hess(4,l)) (Hess(4,2)+Hess(2,4)) (Hess(4,3)+Hess(3,4))

=
(Hess(4,4)) 0 0];
AA(5,:) =
[-1 100 000];
AA(6,:) =
[0001000];
CC(l,:) =
[-dytn(l)];
CC(2,.) =
[-dytn(2)];
CC(3,.) =
[-dytn(3)];
CC(4,:)
=
[-dytn(4)];
CC(5,:) =
[-G(1)];
CC(6,:) =
[-G(3)];
tempmult =
inv(AA)*CC;
dd( 1 1 ) temp_mult( 1,1);

=
,
dd(2,l) =
temp_mult(2,l)
dd(3, 1) =
temp_mult(3,l)
dd(4,l) =
temp_mult(4,l)
u(l) =temp_mult(5,l);
=
159
"multll"
%Filname:

=
(Hess(l,4)+Hess(4,l)) 0]; -10

=
(Hess(2,4)+Hess(4,2)) 0 0 0];

=
(Hess(3,4)+Hess(4,3)) 0 1 0];

=
(Hess(4,4)) 0 0 1];
AA(5,:) =
[-1000000];
AA(6,:) =
[0010000];
AA(7,:) =
[00 0 100 0];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
CC(5,.) =
[-G(1)];
CC(6,:) =
[-G(3)];
CC(7,:) =
[-G(4)];
tempmult =
inv(AA)*CC;
dd( 1 1 )
,
=
temp_mult( 1,1);
dd(2,l) =
temp_mult(2,l);
dd(3,l)
=
temp_mult(3,l);
dd(4, 1 ) =
temp_mult(4, 1 );
160
u(l)=temp_mult(5,l);
=
=
u(4) temp_mult(7, 1);
161
"multl2"
% Filname :

=
(Hess(l,4)+Hess(4,l)) 0]; -1

=
(Hess(2,4)4-Hess(4,2)) 1 0];

=
(Hess(3,4)+Hess(4,3)) 0 0];

=
(Hess(4,4)) 0 0];
AA(5,:) =
[-10000 0];
AA(6,:) =
[010000];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,.) =
[-dytn(3)];
CC(4,:)
=
[-dytn(4)];
CC(5,:) =
[-G(1)];
CC(6,:) =
[-G(2)];
tempmult =
inv(AA)*CC;
dd(l,l) =
temp_mult(l,l);
dd(2,l) =
temp_mult(2,l);
dd(3,l) =
temp_mult(3,l);
dd(4, 1 ) =
temp_mult(4, 1 );
=
162
13"
% Filname : "mult

=
(Hess(l,4)+Hess(4,l)) (Hess(l,5)+Hess(5,l)) 0 0];

-1
AA(2,:) =
[(Hess(2,l)+Hess(l,2)) (Hess(2,2)) (Hess(2,3)+Hess(3,2))
(Hess(2,4)+Hess(4,2)) (Hess(2,5)+Hess(5,2)) 0 1 0];

=

=
(Hess(4,4)) (Hess(4,5)+Hess(5,4)) 0 0 1],
AA(5,:) =
[-1 0000 00];
AA(6,:) =
[0100000];
AA(7,:) =
[0001000];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)],
CC(5,:) =
[-G(1)];
CC(6,.) =
[-G(2)],
CC(7,:)
=
[-G(4)];
tempmult =
inv(AA)*CC;
=
dd(l,l) temp_mult(l,l);
dd(2,l) =
temp_mult(2,l);
=
=
163
=
=
u(4) temp_mult(7, 1 );
164
14"
%Filname: "mult

=
(Hess(l,4)+Hess(4,l)) 0]; -10

=
(Hess(2,4)+Hess(4,2)) 0 1 0];

=
(Hess(3,4)+Hess(4,3)) 0 0 1];

=
(Hess(4,4)) 0 0 0];
AA(6,:) =
[-1000000];
AA(7,:) =
[0100000];
AA(8,:)
=
[00 10 00 0];
CC(l,:) =
[-dytn(l)];
CC(2,:) =
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
CC(5,:) =
[-G(1)];
CC(6,:) =
[-G(2)];
CC(7,:)
=
[-G(3)];
tempmult =
inv(AA)*CC;
dd(l,l) =
temp_mult(l,l);
dd(2,l) =
temp_mult(2,l);
=
dd(4, 1 )
=
temp_mult(4, 1 );
165
=
u(3)=temp_mult(7,l);
166
"multl5"
% Filname :

=
(Hess(l,4)+Hess(4,l)) (Hess(l,5)+Hess(5,l))-100 0];

=
(Hess(2,4)+Hess(4,2)) (Hess(2,5)+Hess<5,2)) 0 10 0];

=

=
(Hess(4,4)) (Hess(4,5)+Hess(5,4)) 0 0 0 1];
AA(5,:) =
[-10000000];
AA(6,:) =
[01000000];
AA(7,:) =
[00010000];
AA(8,:) =
[00001000];
CC(l,:) =
[-dytn(l)];
CC(2,:)
=
[-dytn(2)];
CC(3,:) =
[-dytn(3)];
CC(4,:) =
[-dytn(4)];
CC(5,;) =
[-G(1)];
CC(6,:) =
[-G(2)];
CC(7,:) =
[-G(3)];
CC(8,:)
=
[-G(4)];
tempmult =
inv(AA)*CC;
dd(l,l) =
temp_mult(l,l);
dd(2, 1 ) =
temp_mult(2, 1 );
167
dd(3,l) =
temp_mult(3,l);
dd(4, 1 ) =
temp_mult(4, 1);
=
u(3)=temp_mult(7,l);
=
168
"stepsize"
% filename :
% This subroutine determines the proper step size for the optimization procedure
flag =
0;
Beta = gamma*
(abs_dd)A2;
PHI_o =
Y + R(kk+l)*V;
mm
=
0;
for flag =
0;
t =
(l/2)Amm
X +
=
[xn(l) xn(2) xn(3) xn(4)] t*rranspose(dd)
a(2)
=
X(l)*(a(l)+DELTA);
a(3)
=
X(2)*(a(4)-OH);
Kf =
X(3)*Kf_max;
Kr =
X(4)*Kr_max;
deformationo
PHI_1 =
Y+R(kk+l)*V;
>=
if PHI_1 PHI_o -
t*Beta
=
mm mm+l;
flag =
0;
else
flag=l;
a(2)
=
X(l);
a(3)
=
X(2);
Kf = X(3);
3*9
Kr = X(4);
end
170

Machine Tool Spindle Design PDF

Uploaded by

Copyright:

Available Formats

Machine Tool Spindle Design PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Machine Tool Spindle Design PDF

Uploaded by

Copyright:

Available Formats

Rochester Institute of Technology

RIT Scholar Works

Machine tool spindle design

Follow this and additional works at: http://scholarworks.rit.edu/theses

Dr. Wayne W. Walter

Dr. Hany A. Ghoneim

Dr. Kevin B. Kochersberger

Dr. Satish Kandiikar

DEPARTMENT OF MECHANICAL ENGINEERING

Tide of Thesis: Machine Tool Spindle Design

Date: ~ •/17 /oc Signature of Author: _

In modern machine tool applications the performance of a machine tool is judged

divided into three main sections.

deformed shape of the spindle.

response given the spindle design parameters.

The final section is an optimization of the spindle design. Given constraints on

deflection at the spindle's gauge line.

2.0 Static Analysis 15

2.1 Deflection of Elastic Shaft 17

2.2 Deflection of Bearings 26

2.3 Matlab Solution 27

3.0 Dynamic Analysis 34

3.1 Modal Analysis 35

3.2 Matlab Solution for Mode Shapes 61

3.3 Forced Response 66

3.4 Matlab Solution for Forced Response 70

4.0 Optimization Analysis "4

4.2 The Constrained Steepest Descent Algorithm 81

4.3 The Matlab Solution for Optimization 83

Appendix A : Batch Fiie Template 91

Appendix B: Matiab Programs for Static Analysis 93

Appendix C: Matiab Programs for Dynamic Analysis 107

Appendix D: Matiab Programs for Optimization Analysis 131

3.1 Comparison of Resonant Frequencies (FEA vs. Analytical) 66

4.i Table of Design Variables 76

4.2 Optimum Values for Static Analysis 84

4.3 Optimum Values for Dynamic Analysis 84

2.i Static Spindle Model 16

2.2a Elastic Deflection of Spindle Shaft 18

2.2b Deflection of Spindle Bearings 18

2.3 Model of Uniform Beam 19

2.4 Transformation of Beam Segment 22

2.5 Shear and Bending Moment Diagram for a Beam Segment 23

2.6 Batch File, Static Analysis 28

2.7 Malta Representation of Spindle Geometry 29

2.8a Deflection Contribution of Elastic Shaft 30

2.8b Deflection Contribution of Support Bearings 30

2.9 Total Deflection of Spindle 31

2.10 Comparison of Total Spindle Deflection (FEA vs. Malta) 33

3.1 Dynamic Spindle Model 36

3.2 Lateral Vibration of an Euler-Bernoulli Beam 37

3.3 Sample Euler-Bernoulli Beam 40

3.4 Mode Shapes for Sample Beam 43

3.5 Stepped Euler-Bernoulli Beam 45

3.6 Sample Stepped Euier-Bernouiii Beam 49

3.7a Boundary Conditions for Section i 52

3.7b Boundary Conditions for Section 2 52

3.7d Boundary Conditions for Section 4 53

3.8 Free Body Diagram, Joint 1 55