- Email: [email protected]
A natural solution to the modal equations
Computers Q Structures Vol. 44. No. 6. pp. 1371-l 376. 1992 Printed in Great Britain.
0
0045.7949/92 $5.00 + 0.00 1992 Pergamon Press Ltd
TECHNICAL NOTE A NATURAL
SOLUTION TO THE MODAL EQUATIONS A. E. ANUTA, JR
Civil Engineering Department, Arizona State University, Tempe, AZ 85287, U.S.A. (Receir,ed 19 June 1991) Abstract-The direct solution for all of the mode shapes of undamped, linear elastic systems in free vibration is accomplished by means of a converging iteration scheme applied to the modal equations. Previous versions of the digital computer program used in this solution find successive mode shapes with a decreasing accuracy. The present version employs a natural solution, so that every mode shape is found with the same accuracy.
INTRODUCTION The modal equations are a set of simultaneous cubic equations, obtained for undamped, linear elastic systems in free vibration by minimizing the system energy function. The system energy function is developed through the principle of conservation of energy. The resulting modal coefficients, e, are then given in terms of the system mass and stiffness coefficients, m and k respectively, as e++,= m,k,, - kgm,,
(14
and the corresponding modal equations are [I] i i ~e~k,xjxkx,=O. ,=I*=,,=,
i-l,2
,,...,
n.
equations, wherein every mode shape is found with the same accuracy. THEORY Due to the symmetry of eqn (la), the modal coefficients for repeated indices, such as: e,,,, , e22,2, etc., are always zero, so that each modal equation, i, contains no term involving x:, and this omission facilitates their solution. A relaxation scheme begins with an assumed set of coordinate values, and these starting values are modified incrementally within the problem equations until they approach the solution with an acceptable accuracy. The chosen starting values, x*, are
(lb)
x: = 1, all i
The solution to the modal equations is the complete set of mode shapes, describing all of the free vibration states for the system. The modal equations are always singular, possessing multiple solutions. The sub-set of solutions correctly describing the vibration states is obtained by applying an occult principle: ‘For a singular set of simultaneous equations having a multiplicity of solutions, the only true solution is the occult solution, which occurs when the conflicting requirements of the equations themselves. their spatial and temporal constraints, and any auxiliary equations all fall into a conjunctive alignment’[2]. The auxiliary equations for this problem are the system orthogonality conditions. The solution to the modal equations employs a converging iteration scheme, involving many repetitive calculations, so that it is only practical with the aid of a digital computer. The solution is based upon the Southwell relaxation method [3]; however, the auxiliary equations delimit the paths which the coordinate values are allowed to assume in their progression towards convergence. This results in a ‘forced solution’ to the modal equations. In the existing iteration scheme for solving the modal equations, described in [l], each of the mode shapes is found separately, and it is then used in the orthogonality conditions to force the next mode shape to its convergence. The process contains small compu~tional errors which tend to accumulate, especially in the orthogonahty operation. Thus, the accuracy degrades for each successive mode shape found. The iteration scheme of [l] avoids the issue by allowing the acceptable error at convergence to double for each successive mode shape. It is the purpose of this investigation to develop a ‘natural solution’ to the modal
(28)
and these values, excepting the ith coordinate, are substituted into its corresponding modal equation, which then evolves into a quadratic equation because of the absence of the xf term a,x:2+/3,.v:i-y,=0,
i=I,2
,...,
n,
(2W
where the G(,8, and y coefficients result from the various products of the modal coefficients and the starting coordinate values (i.e., a numerical example illustrating the method is presented in the Appendix of [l]). A changed, or intermediate, value for each coordinate, x**, is obtained by solving eqn (2b) by means of the quadratic formula
x**_-B.*e, I
-
24
i=, tt. 2
,.
,n
(3)
and this intermediate coordinate value immediately replaces the existing value, x *; this operation constitutes one step in the relaxation scheme. The quadratic formula allows a choice of algebraic signs before the radical. If the negative sign is always chosen, the relaxation scheme proceeds in a natural convergence to the system Iowest mode shape, associated with its lowest, or fundamental, frequency. Now, the origin of the modal equations is the system energy function, which for undamped, linear elastic systems in free vibration corresponds to Rayleigh’s quotient [4]. The fundamental mode shape always contains less distortion energy than any of the other mode shapes, so that the exclusive use of the negative sign before the radical in eqn (3) gravitates the solution to
1371
1372
Technical Note
the mode shape containing the least distortion energy in a manner analagous to Rayleigh’s method for estimating the fundamental frequency of a vibrating system [5].t If the positive sign before the radical in eqn (3) is always chosen, the relaxation scheme proceeds in a natural convergence to the system highest mode shape, containing the most distortion energy. The method here is analagous to the Southwell-Dunkerly method for estimating a free vibration frequency [6].$ After the lowest and highest mode shapes are found by the natural solution described above, the remaining n - 2 modes shapes, intermediate to the lowest and highest, are then found successively through an application of the occult principle. The iteration scheme is allowed to choose the algebraic sign before the radical in eqn (3) so as to cause a minimum change in that coordiante value, and the resulting set of intermediate coordinates is orthogonalized against all known mode shapes. The combination of a relaxation step followed by the orthogonalization operation results in a forced solution for the mode shape, with degrading accuracy as it approaches the solution point. From the previous material describing the existing iteration scheme, it is obvious that the method presently employed for the lowest and highest mode shapes is already satisfactory, while the intermediate mode shapes are the cause of the problem under investigation. Also, the successful natural convergences of the lowest and highest mode shapes, purely through the selection of algebraic signs before the radical in eqn (3), suggests that some mixture of positive and negative signs for individual equations might lead to a natural convergence for the intermediate mode shapes. There is a kernel of truth in this idea which is indeed important to the problem solution, but the method by itself is not successful. There are several reasons why a direct application of a mix of sign choices before the radical in eqn (3) fails to always successfutly converge to the correct interm~iat~ mode shape: (1) there is no reliable method for initially identifying the proper set of algebraic signs, (2) a random selection process is impractical because of the large number of sign combinations possible, (3) different solution points can share the same mixture of algebraic signs, and (4) since the intermediate mode shapes represent neither a minimum nor a maximum of distortion energy, as the lowest and highest mode shapes always do, the process often converges to a false solution which is not orthogonal to all the others. Thus, any natural solution for the intermediate mode shapes obviously requires assistance from the orthogonality conditions in order to avoid the false solutions, but then the accumulation of computational errors greatly reduces their effectiveness close to the solution point. A contradiction exists relative to the proper use of the orthogonality conditions in the convergence process. Since the modal equations are singular, the function surface they generate has many minimum points, each of which is a solution point, on a curved function surface. However, it is an established mathematical principle that in the near vicinity of a point, a curved function surface approaches a flat Euclidean surface, From this observation, the philosophy employed in the successful natural solution for the intermediate mode shapes. is to use the existing forced convergence scheme of the occult principle from the beginning of the process down close to a solution point,
t The proof of Rayleigh’s method for estimating the fundamental frequency of an undamped, linear elastic system in free vibration is contained in Rayleigh’s minimum theorem. $ Just as Rayleigh’s method theoretically limits the greatest distortion energy of a vibration frequency, the Southwell-Dunkerly method limits its least distortion energy.
thereby satisfying the orthogonality conditions, and then to abandon it and its accumulated errors for the last few remaining iteration steps to the desired convergence accuracy. The final iteration steps employ only a proper mix of algebraic signs before the radical in eqn (3) in a natural relaxation scheme similar to the one used successfully for the lowest and highest mode shapes. During the iteration cycle just preceding the switch over to the natural convergence scheme, the mixture of algebraic signs which the computer program voluntarily chooses is observed and then held constant afterwards. Close to the solution point, the nearly flat function surface allows the relaxation scheme to successfully converge to a correct solution which is indeed one of the free vibration states of the system. DISCUSSION
A computer program, based upon the original program of [I], modified to employ the natural convergence scheme for ah of the intermediate mode shapes, is presented in the Appendix. There are also some other minor changes incorporated into it which result in further improvements to its performance. Equation (2b) can be put into the form a,x:‘+~,x~+y,=FR6,,
i=l,2
,...,
n,
(4)
where 6 is the equation error term, and FR is an arbitrary fractional number. If FR is zero. eqn (4) is identical to (2b), and if FR is one, d is the existing error for the ith equation. At any solution point for the modal equations, all of the error terms are zero. Thus, the magnitude of the error terms indicates the proximity of the coordinate values to the solution point. ,. After each iteration step, the coordinate values are nor“malized so that the largest among them has an absolute value of one, and in the computer program of the Appendix, a convergence is assumed for any mode shape whenever the first error term reaches an absolute value of less than I x 10W9.Empirical observations suggest that at this definition of convergence, all of the individual coordinates attain a solution accuracy of approximately eight significant figures. This definition of convergence is arbitrary but appears to be consistent with problems caused by the function surface curvature. Solving eqn (4) by means of the quadratic formula gives
,,_-B,~~~i-4a,(i,-FR6,),
x, -
Ia,
i=l
, . 2,
..,n
(5)
and if FR is allowed to have some small negative value, say -0.3 < FR < -0.1, the resulting intermediate coordinate vaiue, x**, exceeds or ‘overshoots’, the value predicted by eqn (3). Usually, this action significantly reduces the number of iteration steps to convergence. The value of FR in this application is arbitrary, and it is optimum at different values for different vibration problems. For the lowest and highest mode shapes, the repeated use of eqn (5) in a relaxation scheme always produces a continuous reduction of the coordinate error terms. The probable reason for this strong performance towards convergence is that these two mode shapes lie on the boundary of the region containing all of the solutions and there is a minimum warping of the function surface here, resulting in a well-behaved convergence. This desirable behavior of the convergence routine is not always observed for the natural solution of the intermediate mode shapes. A possible reason is that the function surface is more irregular in the interior of the region. For whatever reason, it occasionally happens that a change in one of the coordinate values by means of eqn (5) causes the error terms for the other coordinates to increase. It is then possible for the relaxation scheme to stall cyclicly. To avoid this situation, the sum of the absolute
1373
Technical Note values for all n of the error terms is used as the acceptance criterion for any proposed coordinate change. If a change is not acceptable, the program tries several different values of FR, and if none are acceptable, the routine simply moves on to the next coordinate. These repeated trials are time consuming, but fortunately they occur rather infrequently; their use greatly strengthens the convergence process. From the material contained in the Discussion, it is obvious that there are several arbitrary parameters contained in the present version of the computer program: the acceptable error at convergence, ERR, the proximity to the solution where the intermediate mode shapes abandon the forced solution for the natural solution, CRS, and the fractional values used in eqn (S), FR. The values supplied in the present version of the computer program represent conservative values successfully tested in several vibration problems. The singular characteristics of the modal equations are obviously unfavorable to easy problem solutions. Computational difficulties are unavoidable but they are manageable. There are certain advantages to the modal equations method compared to the conventional solution for the vibration problem. Among the advantages: (I) any problem solution can be pursued down to any degree of accuracy deemed desirable, one part in one million is possible provided the computer precision is set to this limit, (2) since
all of the mode shapes are always distinct, the computational difficulties associated with equal or very close frequencies never occur, and (3) because of these two reasons, valid solutions are easily obtained for vibration problems even when the individual mass and/or stiffness ratios are very large; say, equal to or greater than one million to one.
REFERENCES
1. A. Anuta, Jr, The modal equations. Comput. Sfruct. 18, 955-962 (1984). 2. A. Anuta Jr, The mode shape reversal. Camput. Strwt. 33, 106 (1989). 3. R. Southwell, Relaxation Methods in Theoretical Physics, Vol. I, pp. 766. McGraw-Hill, New York (1921). 4. J. Strutt (Lord Rayleigh), The Theory of Sound, 2nd Edn, Vol. I, pp. 1099113. Dover, New York (1945). 5. J. Strutt (Lord Rayleigh), The Theory of Sound, 2nd Edn, Vol. I. pp. 287-289. Dover, New York (1945). 6. R. Southwell, An In~rodue~~o~to rhe Theory of Elasticity, 2nd Edn, Chap. 14. Oxford University Press (1941) and J. Dunkerly, On the whirling and vibration of shafts. Trans R. Sot. AlSJ, 279 (1894).
APPENDIX
The computer program presented here solves the undamped, linear elastic vibration problem for all of the system mode shapes and frequencies. The program is written in FORTRAN IV language. The input vaiues are: the number of coordinates, a program safety-stop number, and the system mass and stiffness matrices. The output includes: the input data, the coordinate starting values, the various mode shapes and their corresponding natural frequencies, the number of iterations required for each mode shape, and the error term at each convergence. c***
c*** C”“R c*** c***
1 2
3
4 5
7” 8 9 c***
MAX IS THE PROGRAM SAFETY STOP, N IS THE NUMBER OF EQUATIONS, A(1, Jf IS THE MASS MATRIX B (I, J) IS THE STIFFNESS MATRIX STOR(I,NS) IS A ROW VECTOR STORAGE FOR ESC(I,J,K,L) NS=(SUM(I*(It1)))/2., FOR I=l,N DIMENSION A(lO,lO),B(lO,lO),~(lO,lO,lO),ES(10,10,1O),~SA(lO,lO),ES 1B(10,10,10),STOR(10,220),ALPKA(10),BETA(10),GAMMA(10),DELTA(10),BO 2T(1O),TOP(lO),G(lO},FREQ(10),SR(lO),X(lO,lO),Y(lO) INTEGER CR. H. PS DOUBLE PRECISION A, B, G, E, ES, ESA,ESB, STOR, ALPHA, BETA, GAMMA,DELTA, SI lGl,SIG2,BOT,TOP,DEN,CRS,ERR,FR,SR,RAD,RN,RP,X,Y READ(S,l)N,MAX READ(5;2) (-~A(I.J),J=l.N).I=l,NI READ(5;2j ((B(I;Jj;J=l;Nj;I=l;Nj FORMAT (215)
FORMAT(5F10.5) DO 3 I=l,N DO 3 J=l,N x(I,J)=l: WRITE(6,4) WRITE (6,5)N,MAX WRITE(6,6) WRITE(6,7) ((A(I,J),J=l,N),I=l,N) WRITE (6,8) WP.ITE(6,7) f(B(I,J),J=l,N),I=l.,N) WRITE (6,s) WRITE (6‘7) ( (X(1, J), J=l,Nf , I=I,N) FORMAT(5X,‘PROGRAM TO FIND ALL OF THE MODE SHAPES AND FREQUENCIES 10F A VIBRATING SYSTEM’//) FORMAT(5X, ‘NUMBER OF EQUATIONS = ’ -14. /.5X. ‘PROGRAM SAFETY STOP AF ITERATIONS’ ,7,5X, ‘THE HIGHEST’ AND LOWEST MODE SHAPES ARE lTER!,I6,‘. 2 FOUND FIRST, THEN THE OTHERS ARE FOUND IN A RANDOM MANNER'//) FORMAT(5X,'M?iSS COEFFICIENTS'/) FORMAT(5X,SE20.9) FOR~T(lH//,SX,'STIFFNESS COEFFICIENTS'/) FO~T(lH//,SX,~STARTING VALUES'/'} THE MODAL EQUATIONS ARE DEVELOPED FROM THE SYSTEM MATRICIES DO 31 I=l,N DO 10 J=l,N DO 10 K=l,N DO 10 L=l,N
1374 10
1': 13 14 15 16 17
Technical Note E(J,K,L)=A(I,J)*B(K,L)-B(I,J)*A(K,L) DO 17 J=l,N DO 17 K=J;N DO 17 L=K.N IF(J-K)12;11,12 IF(K-L)15,14,15 IF(K-L)13,15,13 IF(J-Lb16.15.16 ES(J,K;L)kE(>,K,L) GO TO 17 E; W&&r;) =E (J, K, L) +E (L, J, K) +E (K, L, J)
27
ES(J,K,L)=E(J,K,L)+E(L,J,K)tE(K,L,J)tE(L, CONTINUE PS=O DO 30 J=l,N DO 30 K=J,N DO 30 L=K,N CALL STORE(J,K,L,N,PS,NS) IF(I-J)20,18,20 IF(I-K)19,23,19 IFfI-Lj26.24.26 IFiI-Kj22;21;22 IF(I-L)27,25,27 IF(I-L)29,28,29 ESA(I,L)=ES(J,K,L) GO TO-30 ESA(I,K)=ES(J,K,L) GO TO 30 ESA(I,J)=ES(J,K,L) GO TO 30 ESB(I,K,L)=ES(J,K,L) GO TO 30 ;r$;, $L’=ES (J, K, L)
28
;$W;,$K’=WJ,K,L’
;: z ;z 24 25 26
K, J)tE(K,J,L)tE(J,L,
K)
STOR(I,NS)=ES(J,K,L) CONTINUE CONTINUE z: C*** THE SOLUTION OF THE MODAL EQUATIONS BEGINS ERR=0.000000001 CRS=lO.*ERR DO 83 H=l,N CR=0 M=O 32 M=Mtl IF(MAX-M)71,71,33 IF(II-N)35,34,34 z: II=0 11=11+1 35 I=11 SIGl=O. 36 PS=O ALPHA(I)=O. BETA(I)=O. GAMMA(I)=O. DO 37 J=l,N ALPHA(I)=ALPHA(I)tESA(I,J)*X(H,J) DO 37 K=J,N BETA(I)=BETA(I)tESB(I,J,K)*X(H,J)*X(H,K) DO 37 L=K,N CALL STORE(J,K,L,N,PS,NS) GAMMA(I)=GAMMA(I)+STOR(I,NS)*X(H,J)*X(H,K)*X(H,L) 37 **2tBETA(I)*X(H,I)tGAMMA(I) DELTA(I)=ALPHA(I)*X(H,I) IF(H-2)38,39,41 C*** THE NATURAL SOLUTION OF THE LOWEST AND HIGHEST MODE SHAPES BEGINS 38 SR(I)=l. FR=-0.06 GO TO 40 SR(I)=-1. 39 FR=-0.32 *ALPHA(I)*(GAMMA(I)-FR* 40 X(H,I)=(-BETA(I)+SR(I)*(DSQRT(BETA(I)**2-4. lDELTA(I)))))/(2.*ALPHA(I)) GO TO 67 c*x* THE NATURAL SOLUTION OF THE INTERMEDIATE MODE SHAPES BEGINS DO 42 J=l,N Y(J)=X(H,J) :: IF(M-2*N)58,58,43 IF(CR)44,44,46 IF(DABS(DELTA(II))-CRS)45,58,58 t: CR=1 45 29
Technical Note
51
52 53 54 z"6 57 C*** 58
59 60 C*** 61 62
65 66 CR** 67 68 69
C***
1375
SIGl=SIGltDABS(DELTA(I)) IF(I-N)48,47,47 I=0 1=1+1 IF(I-11)36,49,36 FR=O.lO **2-4.*ALPHA(I)*(GAMMA(I)-FR*DE Y(I)=(-BETA(I)tSR(I)*(DSQRT(BETA(I) lLTA(I)))))/(2.*ALPHA(I)) SIG2=0. PS=O ALPHA(I)=O. BETA(I)=O. GAMMA(I)=O. DO 52 J=l,N ALPHA(I)=ALPHA(I)+ESA(I,J)*Y(J) DO 52 K=J,N BETA(I)=BETA(I)tESB(I,J, K)*Y( J)*Y(K) DO 52 L=K,N CALL STORE(J,K,L,N,PS,NS) GAMMA(I)=GAMMA(I)+STOR(I,NS)*Y(J)*Y(K)*Y(L) SIG2=SIG2+DABS(ALPHA(I)*Y(I)**2+BETA(I)*Y(I)+GAMMA(I)) IF(I-N)S4,53,53 I=0 I=Itl IF(I-11)51,55,51 IF(SIGl-SIG2)56,57,57 FRlFRt0.35 ' IF(FR-1.20)50,32,32 X(H,I)=Y(I) GO TO 67 THE FORCED SOLUTION OF THE INTERMEDIATE MODE SHAPES BEGINS RAD=DSQRT(BETA(I)**2-4. *ALPHA(I)*(GAMMA(I)+0.1S*DELTA(I))) RP=(-BETA(I)tRAD)/(2.*ALPHA(I)) RN=(-BETA(I)-RAD)/(2.*ALPHA(I)) IF(DABS(Y(I)-RP)-DABS(Y(I)-RN))59,S9,60 Y(I)=RP SR(I)=l. GO TO 61 Y(I)=RN SR(I)=-1. THE INTERMEDIATE MODE SHAPES ARE ORTHOGONALIZED DO 62 J=l,N TOP(J)=O. BOT( J)=O. DO 64 L=2,H DO 63 J=l,N DO 63 K=l,N TOP(L)=TOP(L)+X(L-l,J)*A(J,K)*Y (K) BOT(L)=BOT(L)tX(L-l,J)*A(J,K)*X(L-1, K) G(Lj=TOP(Lj/BOTiL) DO 65 J=l,N DO 65 L=2,H Y(J,=Y(J,-G(L)*X(L-l,J) Dd 66 j=i,N X(H,J)=Y(J) ALL THE MODE SHAPES ARE NORMALIZED K=l
DO 69 J=2,N IF((DABS(X(H,K)))-(DABS(X(H,J))))68,69,69 K=J CONTINUE DEN=DABS(X(H,K)) DO 70 J=l,N x(H,J)=x(H,J) /DEN A CONVERGENCE DECISION IS MADE IF(DABS(DELTA(II))-ERR)75,75,32 WRITE(6;72) FORMAT(lH/,SX,'MAXIMUM NUMBER OF ITERATIONS GO TO 75 WRITE(6,74) ~~f$T;;H/,5X,'EQUATION SET IN ERROR')
EXCEEDED')
WRITE(6,76)M,H,DELTA(II) FORMAT(lH//,SX,'THE PROGRAM REQUIRED',I5,' ITERATIONS TO REACH WITH AN ERROR VALUE OF',E17.8,/) 1E SHAPE NUMBER',I3,' WRITB(6,7)(X(H,J),J=l,N) THE FREQUENCIES ARE CALCULATED FROM THE MODE SHAPES AENR=O. BENR=O. DO 77 J=l,N DO 77 K=l,N AENR=AENR+X(H,J)*A(J,K)*X(H,K)
MOD
1376 77
78 79 :: 82 83 C*** 84 a5 86
C""X
100 101 102 103
Technical Note BENR=BENR+X(H,J)*B(J,K)*X(H,K) IF(AENR)73,73,78 TO=4.*AENR*BENR IF(T0)73,79,80 FREQ(H)=O. GO TO 81 F~Q(H)=SQRT(~~/~Z.*~NR) WRITE(6,8Z)H,F~Q(~) FORMAT(lH//,5X,'FREQUENCY NUMBER',I4,' =',E18.8,' RADIANS PER SEC0 lND'//) CONTINUE THE FREQUENCIES ARE PUT INTO AN ASCENDING ORDER J=l DO 86 I=J,N IF(FREQ(J)-FRBQ(I})86,86,85 FRED=FREQ(J) FREQ(J)=FREQ(I) FREQ(I)=FRED CONTINUE J=Jtl IF(J-N)84,84,87 WRITE(6,88) FO~T(~H//,5X,'T~ FREQUENCIES IN ASCENDING ORDER ARE',,') DO 89 I=l.N WRITE(6,8i)I,FREQ(I) WRITE(6,90) FORMAT(lH/,5X,'PROGRAM COMPLETED') STOP END SUBROUTINE STORE IS A ROW VECTOR STORAGE FOR ESC(I,J,K,L) SUBROUTINE STORE(J,K,L,N,PS,NS) INTEGER PS NS=(J-l}*N*Nt(K-l)*N+L-PS IF(L-N)103,100,100 IF(K-L)3.02,101,101 PS=PStJ*(Ntl) GO TO 103 PS=PS+K RETURN END
0
0045.7949/92 $5.00 + 0.00 1992 Pergamon Press Ltd
TECHNICAL NOTE A NATURAL
SOLUTION TO THE MODAL EQUATIONS A. E. ANUTA, JR
Civil Engineering Department, Arizona State University, Tempe, AZ 85287, U.S.A. (Receir,ed 19 June 1991) Abstract-The direct solution for all of the mode shapes of undamped, linear elastic systems in free vibration is accomplished by means of a converging iteration scheme applied to the modal equations. Previous versions of the digital computer program used in this solution find successive mode shapes with a decreasing accuracy. The present version employs a natural solution, so that every mode shape is found with the same accuracy.
INTRODUCTION The modal equations are a set of simultaneous cubic equations, obtained for undamped, linear elastic systems in free vibration by minimizing the system energy function. The system energy function is developed through the principle of conservation of energy. The resulting modal coefficients, e, are then given in terms of the system mass and stiffness coefficients, m and k respectively, as e++,= m,k,, - kgm,,
(14
and the corresponding modal equations are [I] i i ~e~k,xjxkx,=O. ,=I*=,,=,
i-l,2
,,...,
n.
equations, wherein every mode shape is found with the same accuracy. THEORY Due to the symmetry of eqn (la), the modal coefficients for repeated indices, such as: e,,,, , e22,2, etc., are always zero, so that each modal equation, i, contains no term involving x:, and this omission facilitates their solution. A relaxation scheme begins with an assumed set of coordinate values, and these starting values are modified incrementally within the problem equations until they approach the solution with an acceptable accuracy. The chosen starting values, x*, are
(lb)
x: = 1, all i
The solution to the modal equations is the complete set of mode shapes, describing all of the free vibration states for the system. The modal equations are always singular, possessing multiple solutions. The sub-set of solutions correctly describing the vibration states is obtained by applying an occult principle: ‘For a singular set of simultaneous equations having a multiplicity of solutions, the only true solution is the occult solution, which occurs when the conflicting requirements of the equations themselves. their spatial and temporal constraints, and any auxiliary equations all fall into a conjunctive alignment’[2]. The auxiliary equations for this problem are the system orthogonality conditions. The solution to the modal equations employs a converging iteration scheme, involving many repetitive calculations, so that it is only practical with the aid of a digital computer. The solution is based upon the Southwell relaxation method [3]; however, the auxiliary equations delimit the paths which the coordinate values are allowed to assume in their progression towards convergence. This results in a ‘forced solution’ to the modal equations. In the existing iteration scheme for solving the modal equations, described in [l], each of the mode shapes is found separately, and it is then used in the orthogonality conditions to force the next mode shape to its convergence. The process contains small compu~tional errors which tend to accumulate, especially in the orthogonahty operation. Thus, the accuracy degrades for each successive mode shape found. The iteration scheme of [l] avoids the issue by allowing the acceptable error at convergence to double for each successive mode shape. It is the purpose of this investigation to develop a ‘natural solution’ to the modal
(28)
and these values, excepting the ith coordinate, are substituted into its corresponding modal equation, which then evolves into a quadratic equation because of the absence of the xf term a,x:2+/3,.v:i-y,=0,
i=I,2
,...,
n,
(2W
where the G(,8, and y coefficients result from the various products of the modal coefficients and the starting coordinate values (i.e., a numerical example illustrating the method is presented in the Appendix of [l]). A changed, or intermediate, value for each coordinate, x**, is obtained by solving eqn (2b) by means of the quadratic formula
x**_-B.*e, I
-
24
i=, tt. 2
,.
,n
(3)
and this intermediate coordinate value immediately replaces the existing value, x *; this operation constitutes one step in the relaxation scheme. The quadratic formula allows a choice of algebraic signs before the radical. If the negative sign is always chosen, the relaxation scheme proceeds in a natural convergence to the system Iowest mode shape, associated with its lowest, or fundamental, frequency. Now, the origin of the modal equations is the system energy function, which for undamped, linear elastic systems in free vibration corresponds to Rayleigh’s quotient [4]. The fundamental mode shape always contains less distortion energy than any of the other mode shapes, so that the exclusive use of the negative sign before the radical in eqn (3) gravitates the solution to
1371
1372
Technical Note
the mode shape containing the least distortion energy in a manner analagous to Rayleigh’s method for estimating the fundamental frequency of a vibrating system [5].t If the positive sign before the radical in eqn (3) is always chosen, the relaxation scheme proceeds in a natural convergence to the system highest mode shape, containing the most distortion energy. The method here is analagous to the Southwell-Dunkerly method for estimating a free vibration frequency [6].$ After the lowest and highest mode shapes are found by the natural solution described above, the remaining n - 2 modes shapes, intermediate to the lowest and highest, are then found successively through an application of the occult principle. The iteration scheme is allowed to choose the algebraic sign before the radical in eqn (3) so as to cause a minimum change in that coordiante value, and the resulting set of intermediate coordinates is orthogonalized against all known mode shapes. The combination of a relaxation step followed by the orthogonalization operation results in a forced solution for the mode shape, with degrading accuracy as it approaches the solution point. From the previous material describing the existing iteration scheme, it is obvious that the method presently employed for the lowest and highest mode shapes is already satisfactory, while the intermediate mode shapes are the cause of the problem under investigation. Also, the successful natural convergences of the lowest and highest mode shapes, purely through the selection of algebraic signs before the radical in eqn (3), suggests that some mixture of positive and negative signs for individual equations might lead to a natural convergence for the intermediate mode shapes. There is a kernel of truth in this idea which is indeed important to the problem solution, but the method by itself is not successful. There are several reasons why a direct application of a mix of sign choices before the radical in eqn (3) fails to always successfutly converge to the correct interm~iat~ mode shape: (1) there is no reliable method for initially identifying the proper set of algebraic signs, (2) a random selection process is impractical because of the large number of sign combinations possible, (3) different solution points can share the same mixture of algebraic signs, and (4) since the intermediate mode shapes represent neither a minimum nor a maximum of distortion energy, as the lowest and highest mode shapes always do, the process often converges to a false solution which is not orthogonal to all the others. Thus, any natural solution for the intermediate mode shapes obviously requires assistance from the orthogonality conditions in order to avoid the false solutions, but then the accumulation of computational errors greatly reduces their effectiveness close to the solution point. A contradiction exists relative to the proper use of the orthogonality conditions in the convergence process. Since the modal equations are singular, the function surface they generate has many minimum points, each of which is a solution point, on a curved function surface. However, it is an established mathematical principle that in the near vicinity of a point, a curved function surface approaches a flat Euclidean surface, From this observation, the philosophy employed in the successful natural solution for the intermediate mode shapes. is to use the existing forced convergence scheme of the occult principle from the beginning of the process down close to a solution point,
t The proof of Rayleigh’s method for estimating the fundamental frequency of an undamped, linear elastic system in free vibration is contained in Rayleigh’s minimum theorem. $ Just as Rayleigh’s method theoretically limits the greatest distortion energy of a vibration frequency, the Southwell-Dunkerly method limits its least distortion energy.
thereby satisfying the orthogonality conditions, and then to abandon it and its accumulated errors for the last few remaining iteration steps to the desired convergence accuracy. The final iteration steps employ only a proper mix of algebraic signs before the radical in eqn (3) in a natural relaxation scheme similar to the one used successfully for the lowest and highest mode shapes. During the iteration cycle just preceding the switch over to the natural convergence scheme, the mixture of algebraic signs which the computer program voluntarily chooses is observed and then held constant afterwards. Close to the solution point, the nearly flat function surface allows the relaxation scheme to successfully converge to a correct solution which is indeed one of the free vibration states of the system. DISCUSSION
A computer program, based upon the original program of [I], modified to employ the natural convergence scheme for ah of the intermediate mode shapes, is presented in the Appendix. There are also some other minor changes incorporated into it which result in further improvements to its performance. Equation (2b) can be put into the form a,x:‘+~,x~+y,=FR6,,
i=l,2
,...,
n,
(4)
where 6 is the equation error term, and FR is an arbitrary fractional number. If FR is zero. eqn (4) is identical to (2b), and if FR is one, d is the existing error for the ith equation. At any solution point for the modal equations, all of the error terms are zero. Thus, the magnitude of the error terms indicates the proximity of the coordinate values to the solution point. ,. After each iteration step, the coordinate values are nor“malized so that the largest among them has an absolute value of one, and in the computer program of the Appendix, a convergence is assumed for any mode shape whenever the first error term reaches an absolute value of less than I x 10W9.Empirical observations suggest that at this definition of convergence, all of the individual coordinates attain a solution accuracy of approximately eight significant figures. This definition of convergence is arbitrary but appears to be consistent with problems caused by the function surface curvature. Solving eqn (4) by means of the quadratic formula gives
,,_-B,~~~i-4a,(i,-FR6,),
x, -
Ia,
i=l
, . 2,
..,n
(5)
and if FR is allowed to have some small negative value, say -0.3 < FR < -0.1, the resulting intermediate coordinate vaiue, x**, exceeds or ‘overshoots’, the value predicted by eqn (3). Usually, this action significantly reduces the number of iteration steps to convergence. The value of FR in this application is arbitrary, and it is optimum at different values for different vibration problems. For the lowest and highest mode shapes, the repeated use of eqn (5) in a relaxation scheme always produces a continuous reduction of the coordinate error terms. The probable reason for this strong performance towards convergence is that these two mode shapes lie on the boundary of the region containing all of the solutions and there is a minimum warping of the function surface here, resulting in a well-behaved convergence. This desirable behavior of the convergence routine is not always observed for the natural solution of the intermediate mode shapes. A possible reason is that the function surface is more irregular in the interior of the region. For whatever reason, it occasionally happens that a change in one of the coordinate values by means of eqn (5) causes the error terms for the other coordinates to increase. It is then possible for the relaxation scheme to stall cyclicly. To avoid this situation, the sum of the absolute
1373
Technical Note values for all n of the error terms is used as the acceptance criterion for any proposed coordinate change. If a change is not acceptable, the program tries several different values of FR, and if none are acceptable, the routine simply moves on to the next coordinate. These repeated trials are time consuming, but fortunately they occur rather infrequently; their use greatly strengthens the convergence process. From the material contained in the Discussion, it is obvious that there are several arbitrary parameters contained in the present version of the computer program: the acceptable error at convergence, ERR, the proximity to the solution where the intermediate mode shapes abandon the forced solution for the natural solution, CRS, and the fractional values used in eqn (S), FR. The values supplied in the present version of the computer program represent conservative values successfully tested in several vibration problems. The singular characteristics of the modal equations are obviously unfavorable to easy problem solutions. Computational difficulties are unavoidable but they are manageable. There are certain advantages to the modal equations method compared to the conventional solution for the vibration problem. Among the advantages: (I) any problem solution can be pursued down to any degree of accuracy deemed desirable, one part in one million is possible provided the computer precision is set to this limit, (2) since
all of the mode shapes are always distinct, the computational difficulties associated with equal or very close frequencies never occur, and (3) because of these two reasons, valid solutions are easily obtained for vibration problems even when the individual mass and/or stiffness ratios are very large; say, equal to or greater than one million to one.
REFERENCES
1. A. Anuta, Jr, The modal equations. Comput. Sfruct. 18, 955-962 (1984). 2. A. Anuta Jr, The mode shape reversal. Camput. Strwt. 33, 106 (1989). 3. R. Southwell, Relaxation Methods in Theoretical Physics, Vol. I, pp. 766. McGraw-Hill, New York (1921). 4. J. Strutt (Lord Rayleigh), The Theory of Sound, 2nd Edn, Vol. I, pp. 1099113. Dover, New York (1945). 5. J. Strutt (Lord Rayleigh), The Theory of Sound, 2nd Edn, Vol. I. pp. 287-289. Dover, New York (1945). 6. R. Southwell, An In~rodue~~o~to rhe Theory of Elasticity, 2nd Edn, Chap. 14. Oxford University Press (1941) and J. Dunkerly, On the whirling and vibration of shafts. Trans R. Sot. AlSJ, 279 (1894).
APPENDIX
The computer program presented here solves the undamped, linear elastic vibration problem for all of the system mode shapes and frequencies. The program is written in FORTRAN IV language. The input vaiues are: the number of coordinates, a program safety-stop number, and the system mass and stiffness matrices. The output includes: the input data, the coordinate starting values, the various mode shapes and their corresponding natural frequencies, the number of iterations required for each mode shape, and the error term at each convergence. c***
c*** C”“R c*** c***
1 2
3
4 5
7” 8 9 c***
MAX IS THE PROGRAM SAFETY STOP, N IS THE NUMBER OF EQUATIONS, A(1, Jf IS THE MASS MATRIX B (I, J) IS THE STIFFNESS MATRIX STOR(I,NS) IS A ROW VECTOR STORAGE FOR ESC(I,J,K,L) NS=(SUM(I*(It1)))/2., FOR I=l,N DIMENSION A(lO,lO),B(lO,lO),~(lO,lO,lO),ES(10,10,1O),~SA(lO,lO),ES 1B(10,10,10),STOR(10,220),ALPKA(10),BETA(10),GAMMA(10),DELTA(10),BO 2T(1O),TOP(lO),G(lO},FREQ(10),SR(lO),X(lO,lO),Y(lO) INTEGER CR. H. PS DOUBLE PRECISION A, B, G, E, ES, ESA,ESB, STOR, ALPHA, BETA, GAMMA,DELTA, SI lGl,SIG2,BOT,TOP,DEN,CRS,ERR,FR,SR,RAD,RN,RP,X,Y READ(S,l)N,MAX READ(5;2) (-~A(I.J),J=l.N).I=l,NI READ(5;2j ((B(I;Jj;J=l;Nj;I=l;Nj FORMAT (215)
FORMAT(5F10.5) DO 3 I=l,N DO 3 J=l,N x(I,J)=l: WRITE(6,4) WRITE (6,5)N,MAX WRITE(6,6) WRITE(6,7) ((A(I,J),J=l,N),I=l,N) WRITE (6,8) WP.ITE(6,7) f(B(I,J),J=l,N),I=l.,N) WRITE (6,s) WRITE (6‘7) ( (X(1, J), J=l,Nf , I=I,N) FORMAT(5X,‘PROGRAM TO FIND ALL OF THE MODE SHAPES AND FREQUENCIES 10F A VIBRATING SYSTEM’//) FORMAT(5X, ‘NUMBER OF EQUATIONS = ’ -14. /.5X. ‘PROGRAM SAFETY STOP AF ITERATIONS’ ,7,5X, ‘THE HIGHEST’ AND LOWEST MODE SHAPES ARE lTER!,I6,‘. 2 FOUND FIRST, THEN THE OTHERS ARE FOUND IN A RANDOM MANNER'//) FORMAT(5X,'M?iSS COEFFICIENTS'/) FORMAT(5X,SE20.9) FOR~T(lH//,SX,'STIFFNESS COEFFICIENTS'/) FO~T(lH//,SX,~STARTING VALUES'/'} THE MODAL EQUATIONS ARE DEVELOPED FROM THE SYSTEM MATRICIES DO 31 I=l,N DO 10 J=l,N DO 10 K=l,N DO 10 L=l,N
1374 10
1': 13 14 15 16 17
Technical Note E(J,K,L)=A(I,J)*B(K,L)-B(I,J)*A(K,L) DO 17 J=l,N DO 17 K=J;N DO 17 L=K.N IF(J-K)12;11,12 IF(K-L)15,14,15 IF(K-L)13,15,13 IF(J-Lb16.15.16 ES(J,K;L)kE(>,K,L) GO TO 17 E; W&&r;) =E (J, K, L) +E (L, J, K) +E (K, L, J)
27
ES(J,K,L)=E(J,K,L)+E(L,J,K)tE(K,L,J)tE(L, CONTINUE PS=O DO 30 J=l,N DO 30 K=J,N DO 30 L=K,N CALL STORE(J,K,L,N,PS,NS) IF(I-J)20,18,20 IF(I-K)19,23,19 IFfI-Lj26.24.26 IFiI-Kj22;21;22 IF(I-L)27,25,27 IF(I-L)29,28,29 ESA(I,L)=ES(J,K,L) GO TO-30 ESA(I,K)=ES(J,K,L) GO TO 30 ESA(I,J)=ES(J,K,L) GO TO 30 ESB(I,K,L)=ES(J,K,L) GO TO 30 ;r$;, $L’=ES (J, K, L)
28
;$W;,$K’=WJ,K,L’
;: z ;z 24 25 26
K, J)tE(K,J,L)tE(J,L,
K)
STOR(I,NS)=ES(J,K,L) CONTINUE CONTINUE z: C*** THE SOLUTION OF THE MODAL EQUATIONS BEGINS ERR=0.000000001 CRS=lO.*ERR DO 83 H=l,N CR=0 M=O 32 M=Mtl IF(MAX-M)71,71,33 IF(II-N)35,34,34 z: II=0 11=11+1 35 I=11 SIGl=O. 36 PS=O ALPHA(I)=O. BETA(I)=O. GAMMA(I)=O. DO 37 J=l,N ALPHA(I)=ALPHA(I)tESA(I,J)*X(H,J) DO 37 K=J,N BETA(I)=BETA(I)tESB(I,J,K)*X(H,J)*X(H,K) DO 37 L=K,N CALL STORE(J,K,L,N,PS,NS) GAMMA(I)=GAMMA(I)+STOR(I,NS)*X(H,J)*X(H,K)*X(H,L) 37 **2tBETA(I)*X(H,I)tGAMMA(I) DELTA(I)=ALPHA(I)*X(H,I) IF(H-2)38,39,41 C*** THE NATURAL SOLUTION OF THE LOWEST AND HIGHEST MODE SHAPES BEGINS 38 SR(I)=l. FR=-0.06 GO TO 40 SR(I)=-1. 39 FR=-0.32 *ALPHA(I)*(GAMMA(I)-FR* 40 X(H,I)=(-BETA(I)+SR(I)*(DSQRT(BETA(I)**2-4. lDELTA(I)))))/(2.*ALPHA(I)) GO TO 67 c*x* THE NATURAL SOLUTION OF THE INTERMEDIATE MODE SHAPES BEGINS DO 42 J=l,N Y(J)=X(H,J) :: IF(M-2*N)58,58,43 IF(CR)44,44,46 IF(DABS(DELTA(II))-CRS)45,58,58 t: CR=1 45 29
Technical Note
51
52 53 54 z"6 57 C*** 58
59 60 C*** 61 62
65 66 CR** 67 68 69
C***
1375
SIGl=SIGltDABS(DELTA(I)) IF(I-N)48,47,47 I=0 1=1+1 IF(I-11)36,49,36 FR=O.lO **2-4.*ALPHA(I)*(GAMMA(I)-FR*DE Y(I)=(-BETA(I)tSR(I)*(DSQRT(BETA(I) lLTA(I)))))/(2.*ALPHA(I)) SIG2=0. PS=O ALPHA(I)=O. BETA(I)=O. GAMMA(I)=O. DO 52 J=l,N ALPHA(I)=ALPHA(I)+ESA(I,J)*Y(J) DO 52 K=J,N BETA(I)=BETA(I)tESB(I,J, K)*Y( J)*Y(K) DO 52 L=K,N CALL STORE(J,K,L,N,PS,NS) GAMMA(I)=GAMMA(I)+STOR(I,NS)*Y(J)*Y(K)*Y(L) SIG2=SIG2+DABS(ALPHA(I)*Y(I)**2+BETA(I)*Y(I)+GAMMA(I)) IF(I-N)S4,53,53 I=0 I=Itl IF(I-11)51,55,51 IF(SIGl-SIG2)56,57,57 FRlFRt0.35 ' IF(FR-1.20)50,32,32 X(H,I)=Y(I) GO TO 67 THE FORCED SOLUTION OF THE INTERMEDIATE MODE SHAPES BEGINS RAD=DSQRT(BETA(I)**2-4. *ALPHA(I)*(GAMMA(I)+0.1S*DELTA(I))) RP=(-BETA(I)tRAD)/(2.*ALPHA(I)) RN=(-BETA(I)-RAD)/(2.*ALPHA(I)) IF(DABS(Y(I)-RP)-DABS(Y(I)-RN))59,S9,60 Y(I)=RP SR(I)=l. GO TO 61 Y(I)=RN SR(I)=-1. THE INTERMEDIATE MODE SHAPES ARE ORTHOGONALIZED DO 62 J=l,N TOP(J)=O. BOT( J)=O. DO 64 L=2,H DO 63 J=l,N DO 63 K=l,N TOP(L)=TOP(L)+X(L-l,J)*A(J,K)*Y (K) BOT(L)=BOT(L)tX(L-l,J)*A(J,K)*X(L-1, K) G(Lj=TOP(Lj/BOTiL) DO 65 J=l,N DO 65 L=2,H Y(J,=Y(J,-G(L)*X(L-l,J) Dd 66 j=i,N X(H,J)=Y(J) ALL THE MODE SHAPES ARE NORMALIZED K=l
DO 69 J=2,N IF((DABS(X(H,K)))-(DABS(X(H,J))))68,69,69 K=J CONTINUE DEN=DABS(X(H,K)) DO 70 J=l,N x(H,J)=x(H,J) /DEN A CONVERGENCE DECISION IS MADE IF(DABS(DELTA(II))-ERR)75,75,32 WRITE(6;72) FORMAT(lH/,SX,'MAXIMUM NUMBER OF ITERATIONS GO TO 75 WRITE(6,74) ~~f$T;;H/,5X,'EQUATION SET IN ERROR')
EXCEEDED')
WRITE(6,76)M,H,DELTA(II) FORMAT(lH//,SX,'THE PROGRAM REQUIRED',I5,' ITERATIONS TO REACH WITH AN ERROR VALUE OF',E17.8,/) 1E SHAPE NUMBER',I3,' WRITB(6,7)(X(H,J),J=l,N) THE FREQUENCIES ARE CALCULATED FROM THE MODE SHAPES AENR=O. BENR=O. DO 77 J=l,N DO 77 K=l,N AENR=AENR+X(H,J)*A(J,K)*X(H,K)
MOD
1376 77
78 79 :: 82 83 C*** 84 a5 86
C""X
100 101 102 103
Technical Note BENR=BENR+X(H,J)*B(J,K)*X(H,K) IF(AENR)73,73,78 TO=4.*AENR*BENR IF(T0)73,79,80 FREQ(H)=O. GO TO 81 F~Q(H)=SQRT(~~/~Z.*~NR) WRITE(6,8Z)H,F~Q(~) FORMAT(lH//,5X,'FREQUENCY NUMBER',I4,' =',E18.8,' RADIANS PER SEC0 lND'//) CONTINUE THE FREQUENCIES ARE PUT INTO AN ASCENDING ORDER J=l DO 86 I=J,N IF(FREQ(J)-FRBQ(I})86,86,85 FRED=FREQ(J) FREQ(J)=FREQ(I) FREQ(I)=FRED CONTINUE J=Jtl IF(J-N)84,84,87 WRITE(6,88) FO~T(~H//,5X,'T~ FREQUENCIES IN ASCENDING ORDER ARE',,') DO 89 I=l.N WRITE(6,8i)I,FREQ(I) WRITE(6,90) FORMAT(lH/,5X,'PROGRAM COMPLETED') STOP END SUBROUTINE STORE IS A ROW VECTOR STORAGE FOR ESC(I,J,K,L) SUBROUTINE STORE(J,K,L,N,PS,NS) INTEGER PS NS=(J-l}*N*Nt(K-l)*N+L-PS IF(L-N)103,100,100 IF(K-L)3.02,101,101 PS=PStJ*(Ntl) GO TO 103 PS=PS+K RETURN END