C3●

LAMS-2954

.

elC-14 REPORT COLLECTION

REPRODUCTION

COPY

LOS ALAMOS SCIENTIFIC LABORATORY

OF THEUNIVERSITYOF CALIFORNIAo LOSALAMOS NEW MEXICO

COMPUTING THE CRITICAL FREQUENCIES

OF STEPPED SHAFTS

,_ –— .

<=

—=

-

1

!!$-------- —...-.

iCTJ __

-—— ———
LEGAL NOTICE

This reportwas preparedas an accountof Government t3ponsoredwork. Neitherthe UnitedStates,nor the

Commission,norany personactingon behalfof the Commission:

A. Makes anywarrantyor representation, expressed

or implied,withrespecttotheaccuracy,completeness,or

usefulnessof the informationcontainedin thisreport,or

thattheuse of any information,apparatus,method,or process disclosedin thisreportmay notinfringeprivately

owned rights;or .

B. Assumes any liabilities withrespectto theuse

of,or fordamages resultingfrom theuse of any information,apparatus,method,or process disclosedin thisreport.

As used in theabove,“personactingon behalfof the

Commission” includesany employeeor contractorof the

Commission,or employeeof suchcontractor,to theextent

thatsuch employee or contractorof the Commission,or

employee of such contractorprepares,disseminates,or

providesaccessto,any informationpursuantto his employmentor contractwiththe Commission,or hisemployment withsuch contractor.

Printedin USA. Price $0.50. Availablefrom the

Officeof TechnicalServices

U. S. Departmentof Commerce

Washington25,D. C.
LAMs-2954

uC-32, MATHEMATICS AND

COMPUTERS

TID-4500 (22ndEd.)

LOS ALAMOS SCIENTIFIC LABORATORY

OF THE UNIVERSITYOF CALIFORNIA LOS ALAMOS NEW MEXICO

REPORT WRITTEN: August 19, 1963

REPORT DISTRIBUTED: October 11, 1963

COMPUTING THE CRITICAL FREQUENCIES

OF STEPPED SHAFTS

by

BurtonWendroff

=9= M--’- ,r——

,.

, II

-m ~ Contract —— V

W-7405-ENG. 36 withthe u. s. Atomic Energy Commission

All LAMS reperte are informal deeumeuts, ueufdly prepared for a special purpose and primarily prepared for use wftbln the Laboratory rather than for

genexal distribution. This repert baa net been edited, reviewed, or verified

for accuracy. All LAMS reports express the vIewe of tba authora aa of the

time they were wrftten and do net naceeswlly refleet tie opfxdonaof tbe Los

Alamos Scientific Laboratory or tbe final opinion of the authors en the subject.

1
ABSTRACT

Two finite difference methods are presented for representing the

t

critxicsl frequencies of a shaft as the eigenvalues of a matrix. The

matrices are well suited for high speed digitsl. compu-tation. The

methods are applied to a uniform shaft end compared with the known

frequencies.

ACKNOWLEDGEMENT

We would like to thank L. Linson for his help on this problem.
I. A simple model of the steady displacement of a rotating shaft

is provided by the differential equation

(1) (EIy”)“ = A2qy, EI>O, q>0

with boundary conditions; for example, if the shaft is simply supported,

y(0) = y“(0) = y(L) = y“(L) where L is the length of the shaft. ,EI is

the flexural rigidity and q is the linear mass density. The eigenvdues

of this differential problem are known as the critical frequencies

the shaft.

The literature on the subject of numerical approximation of X2

seems to be devoted to the adaptation of Mikelstadt?s method to high

speed computing, with one exception, that being a recent paper [2] on

an integral equation approach. Mikelstadt’s method is described in [3].

Neither of these methods uses a direct attack on (l). In this paper we

apply two well-known and powerful finite difference methods directly to

(1), obtaining in both cases linear systems of the form AX = V2BX Were

A and B are symmetfic, positive definite, and striped. In the first case

B is diagonal and ~~ are obtained by the LR method. In the second case

we evaluate detlA - V2BI by Gaussian elimination and find its zeroes by

3
regula falsi. The first method provides apparent lower bounds for A2;

the second gives rigorous u~er bounds.

II. The first method is based on a continuity principle. To

simplify the notation let p = EI, so that

(2) (PY”)” = ~2u

We suppose Mat p and q are step functions with

continuity. If we now assume that solutions of

then by successive integrations it can be shown

functions are continuous

hl= (PY”)’

ha

= PY”

he = y,

h4=y

common points of dis-

(2) are integrable,

that the following

It must be understood that these conttiuity conditions are physicsl

ones. If they do not hold it means that (1) is not the proper model.

4
We set up a mesh

‘%= ’$C-l +&/, ‘= 0$ 1) ““”> N

x= o

0

in such a way that every discontinuity of p and q is some ~. Set

‘k+1/2 = P(~=/2 ) , etc.

The derivation of the difference equations is started by setting

%-1/2

2

yk =

\

(PY”)” ~

%+=/2+ %./2%1,2 q

‘@Grq22+h’:.i’)

If we approximate the continuity of y by the nelation

‘k

I

%1/2

(PY”)” ~= 1

I

titi

%;/2 q

%+1/2 q

%1/2 %

5
the. we can elimi.tute ~ to obt.i.

where

qk+l/= %1-l/2+ %42%1/2

>= 2

Since all functions are amooti at k+l/2, put

Proceeding as before, we write

%1/2

%= 2

%1/2 + %+1/2 /

py” dx

%1/2

%1/2

6

‘i
to get

where

2Pk-z/a ‘lstl/2

~k =

‘k-l/z %+1/2 + ‘k+l/2 %-1/2

Finally

‘k+l - Yk

%1/2 = ++1,2

The resultant difference equation is

[(

3k!t2 - Y~L ‘k+l

)

- yk

(2a) ~2 6kyk =

%:1/2

‘k+ 1

%3/2 - %+1/2

(

%1

)]

- Yk Yk - Yk-l

- ~k

%+1/2 - 442

[J%1 )

-

Yk Yk - Yk-=

-—

+:1/2

i

%1/2 - %-1/2

.

(

Yk - Yk-1

- zk-l

442 -

‘k-1

)

-‘k-2

442

7
with boundary conditions

Y. =yN. o

P. = %=0

In matrix form we have

(3) Ax= &2BX

where

A= {a],ij B = [bi3}

ai j = aji

ai j =Oi+-jl >2

1

aii = A

i+ 1/2 (

1-1+

=i+~ ~ —

i+l/2

‘pi A

i+1/2

.

1

‘i- 1/2)
‘Aj.~/~[Fi(-+++ =i+l-]

ai’ i-i-l= ‘-[5i+l.(Ai~~/* +-)+5i-]

1

=, A

1 —.

+ A.

L- 1/2 i+1/2

1 1

ai’ i+2 = A ‘i+l Ai+=12 i+l/2

bij=o~ q=s

b

ii = (<)

III. The second method is based on a variation&1 principle. Let

L

[

P(Y”)2 M

(4) -6

P2(Y) =~

Let Z be the space of &U. continuously differentiable functions on

L with piecewise continuous second derivatives which slso satisfy the

9
boundary conditions. Let Smbe any m-dimension&1 subspace

~, i = 1, 2, ‘CO, m be the stationary vslue of V2(y) as y

J.

sm. Then it is known that

X2 <ri, i= 1, 2, ●.”, m

iWe choose as the subspace C2( N-1)$ which we

all continuous functions having continuous first

are c~ic polynomials in the intervsls (~ , %1

define as

of z. Let

varies wer

the set of

derivatives, and which

),k=O, 1,-~*, N-l.

‘ach ‘iCtion h C2 (N- 1) is determined by the value of it and its first

derivative at all the interior mesh points. CsJJ these values yi, and

Y& i =1, 2, ●**, N- 1. Then for any y(x) in C2(N =,

y(x)

‘Theconditions

=Y(x> Y~$ ““0 ) Y~> ““”)

for stationarity of V2(y) are

w Y) ~=~’ ‘=1’ 2’”””’ N-1

or

10
and

L

P L

L

so that 8gain

(6) Ax = C2BX

u.
where A and B are of necessity positive definite snd can be chosen

symmetric. The calculation of A and B is extremely tedious and will

be omitted. We will.give the equations for the cubics in order to aid

the reader who wishes to use this method with other boundary conditions.

We have

y(x) = ax3+b#+cx+d

where

a) For O<i<i+l<N $xi _< x ~ ‘i+l

= Ai+l/2(q + Y;+J - 2(Yf+x - Yi)

a

@i+l/J3

3(Yi+l - Y+ - A~+=/Jq + Y$+.J

b =

(Ai+l/2)2

b) In the simply supported case, for O < x < xl --

12
bO=

3Yi - Y~ Al/2

c =

2A1/2

dO=

3yN-= + +-1/2 ‘h

c =-

%-1/2

b=d=O

If the vector x is taken to be

13

—
then the matrices A and B are block tridiagonsl, as follows:

(\\

b= Cl

a2 b2 c= o\

B =

o

%?-1bN-l/

I ‘1 ‘24Pi-=/, ‘l*i 3./2

A?

1- 1/2 ‘?- 1/2

a. =

1

Qi-~/~ 4pi-1/2

‘:- 1/2 A.1-1/2
/

24pi+l/2

‘:+ 1/2

[

A?

1+ 1/2

‘l~i+l/2

\

%+ 1/2

(

24PSJ2

%/2

%/2

4/2

A=

i- 1/2

12pi-112

A?

1- 1/2

\

l@i-./2

% 1/2

1

/

4pi+l/2

‘i+ 1/2

+ 6%/2

A2

1/2

A?1+1/2

A?

1- 1/2

8pi+~/2 8pi-1/2

A

i+1/2

/

‘i- 1/2

%/2 +

A

3/2 )

%/2

A2

1/2

%!S

A

/

1/2

15
\%,

2

-3 2 %1/2

/

8PN42 + 6pIV.1/2

%3/2 %-1/2

16
Iv.To solve (3) for~= we use the SST method [1]. First, since

B is diagonal with positive diagonsl elements we can write

17
me matrices A and B are positive definite and symmetric, but B-lA

is not symmetric. It can be symmetrized by setting

so that B-l= ~2

Then if

D = CAC

we have

CDC-l = C2A = B-lA

so that D is similar to B-LA and has the ssme eigenvalues. D iS

symmetric and positive definite. Any such matrix D has a square root;

that is,there exists a matrix S, lower triangular, such that

18

D = SST
If D is striped so is S. The SST algorithm is the following:

a) Let D(0)=D

b) Given D(k), find S(k) s~~hthat

~(k) = S(k)S(k)T

c) Let D(k+l)=S(k)TS(k)

The sequence of matrices D

(k)

converges to a diagonal matrix with the

eigenvslue of D(0) onthediagonde ‘l%eoff-diagon ~element so fthe

last row and column converge to zero first, at which time they C= be

discarded and the process continued with successively smaller matrices.

For our particular .(0) the smsllest eigenvalue appears first, then

the next

The

smsllest, etc.

matrix S is found by the following recursion. If D= (Dij)

r

j-l

s. =

Jj ‘Js - x

S2

jr

1=1

(

j-1

E )()-1

‘ij = ‘ij -

s

ir ‘jr ‘jj

*1

i= J+l, j+2,000

s =1, 2, ““”

19
If ~ now consider the variation~ method ~ find that the ~trices

occurring there do not lend themselves to the above method, the reason

being that now B is not diagonsl so B-lA is not striped. It turns out

to be quite feasible with the Stretch computer to find the zeroes of

detlA - TBI, with A and B defined in part III, by direct computation.

So, let

f(T) = detlA - ‘rB1

For any T

pivoting,

the determinant is evsluated by Gaussian elimination without

but programmed to take advantage of the sparseness of A - ‘rB.

The zeroes of f(-r)

~(k+l) =

are found by regula falsi; that is, the iteration

T(k)- f(~k)

~k - ~k-1

f(Tk) - f(’rk-l)

To find ~i we start with

o

T = g

r’ V:(1+ c)

for some smsll 6. The ~~ are the solutions of (3)0

We are able to report on only one actusl calculation with the

20
ahove method. We took

P=q=l

L = 10

We used 20 intervsls with the continuity method and 10 intervsls with

the variational method

~02 :

and obtained the following values:

1.56757 1.570799 1.57081

6.23168 6.2832 6.28387

13.877 14.137 14.145

For this problem it is easy to show that ~ < A; however, this is

not known in genersl. We have observed that for some problems with

discontinuous p and q the w?s seem to increase with decreasing Ax.

St is clear from the table that the upper bounds ~ are more

accurate than the ~; however, this comparison is not completely fair

since the matrices determining ~ have seven non-zero diagonals, while

those detemnining p have five. The

coxupsredwith a seven point form of

not been able to work this form out

variational method should be

the continuity method, but we have

in the genersllcase.

21
v. We should like to conclude with a proof that A in (3) is

positive definite.

Forsnyvectorw= (yl, ● 00 , yN-l) we have to show that (Aw,w) >0

if w # 0, where

N-1

Define the following operators

(Eu)i= Ui+=,2

(k)i = (E- E-l)U= Ui+l/Z - ‘ii/2

The operator E is unitary on the set of UTS which vanish outside

some titervsl, which we may represent by ~ = O if k ~ O or k ~ N. The

unitarity means

& = E-L , where E* = adjoint of E

!lMerefore

(:U,V) = .(U,2V)

22
Let

h

1

i+L/2 = A.

1+1/2

An examination of (2a) shows that

AA AA

Aw = (&&~)W

Putting R = z

we have

R4t=R

Therefore

(Aw,w) = (R@w,w) = (~Rw,Rw)

Equslity can occur only if Rw = O, which means ~w is constant; however,

since W. = WN = 0, this implies w = O.

23
Bibliography

10 Rutishauser, H., “Solution of Eigenvalue Problems with LR transformation,” in NBS Appl. Math. Series, No. 49, 1958, p. k7.

2. Stern, M. and Stippes, M., A unified numerical treatment for

vibrating beams, J. Sot. Ind. Appl. Math, v. U., No. 2, 413 (1963).

3. Timoshenko, S., and D. H. Young, Vibration Problems in Engineering,

Third Ed., D. Van Nostrand Company, Inc. 1955.

24