An efficient approach for whirling speeds and mode shapes of uniform and nonuniform Timoshenko shafts mounted by arbitrary rigid disks

In theory, the whirling motion of a shaft‐disk system is three‐dimensional (3D), however, if the transverse displacement in the vertical principal xy‐plane and that in the horizontal principal xz‐plane for the cross‐section of the shaft located at the axial coordinate x are represented by a complex number, then the mass moment of inertia for unit shaft length (at x) and that for each rigid disk i can be combined with their associated gyroscopic moments (GMs) to form the frequency‐dependent equivalent mass moments of inertia, respectively, in the equations of motion for the rotating Timoshenko shaft carrying arbitrary rigid disks. It is found that the above‐mentioned equations for the rotating 3D shaft‐disk system take the same forms as the corresponding ones for the associated stationary two‐dimensional (2D) Timoshenko beam carrying the same number of disks, so that the approaches for the free vibration analyses of the stationary 2D beam‐disk system can be used to solve the title problem for obtaining the whirling speeds and mode shapes of a rotating 3D shaft‐disk system. Since the order of the eigenproblem equation derived from the presented approach is much smaller than that derived from the conventional finite element method (FEM), the computer time consumed by the former is much less than that by the latter. In addition, the solutions obtained from the presented approach are exact and may be the benchmark for evaluating the accuracy of solutions of the other approximate methods. Numerical examples reveal that the presented approach is available for the uniform or nonuniform shaft‐disk system, and the obtained results are very close to those obtained from existing literature or the FEM. The formulation of this article is available for a shaft‐disk system with various boundary conditions (BCs), to save space, only the cases with pinned‐pinned BCs are illustrated.


INTRODUCTION
The critical speeds of flexible rotors are the important information for the engineers concerned, thus a lot of researchers devoted themselves to the studies of this topic. For example, Prohl 1 calculated the critical speeds of flexible rotors by using a step-by-step integration process analogous to method of Holzer 2 for obtaining the natural frequencies of a torsional system. Later, Myklestad 3 solved the similar problem by using the transfer matrix method (TMM). Neglecting the shaft mass, Green 4 studied the gyroscopic effect of disks on the critical speed of a shaft-disk system. To improve the accuracy of the existing methods, Eshleman and Eubanks 5 presented an analytical method to obtain the first four forward and backward whirling speeds of a uniform Euler-Bernoulli shaft (or Euler shaft, for simplicity) carrying a single rigid disk with the distributed mass of the shaft considered. Two years later, Eshleman and Eubanks 6 studied the critical speeds of a uniform bare Timoshenko shaft (without disks) by using the analytical method with effects of shear deformation (SD), rotary inertia (RI), and gyroscopic moment (GM) of the shaft considered. Shiau and Hwang 7 employed the assumed mode method to study the critical speeds of a rotor-bearing system. By neglecting the gyroscopic effects, Firoozian and Zhu 8 presented a hybrid method for the vibration analysis of the rotor-bearing systems. Zu and Han 9 determined the natural frequencies and associated mode shapes of a spin uniform bare Timoshenko beam with general boundary conditions (BCs) by using the conventional analytical method. Han and Zu 10 extended their previous study 9 to a spin Timoshenko beam subjected to a moving load. With the rotor modeled by a series of distributed and lumped elements, Aleyaasin et al 11 studied the flexural vibration of rotors mounted on fluid film bearings by using the frequency-dependent TMM. Wu et al 12 obtained the exact solution for the whirling speeds and mode shapes of a distributed-mass shaft carrying arbitrary rigid disks with the analytical method. Firouzi and Ghassemi 13 studied the effects of diameter and torsional stiffness on the whirling speed of the ship shafting system with GM of the rotors neglected by analytical method. Bharti et al 14 studied the Sommerfeld effect (the jump phenomenon) at forward and backward critical speeds in a (four degree of freedom) "rigid" rotor shaft system with anisotropic supports by using the analytical method. In addition to the above-mentioned analytical methods or TMMs, Nelson and McVaugh 15 studied the dynamics of the rotor-bearing systems with the finite element method (FEM). Since the last FEM does not consider the effects of SD and RI of the shaft elements, Nelson 16 further presented a finite rotating shaft element using the Timoshenko beam theory. Yamamoto and Ishida 17 introduced the applications of analytical methods, TMM and FEM to the linear and nonlinear dynamics of multidisk rotor-bearing systems. Irani Rahaghi et al 18 studied the longitudinal-torsional and two plane transverse vibrations of a composite Timoshenko rotor by using the differential quadrature method, and two kinds of critical speeds are considered for the rotor: the resonance speed and the instability speed. Torabi and Afshari 19 obtained the exact solution for whirling analysis of axial-loaded Timoshenko rotor by using the basic functions, where the effect of spin speed, axial load, slenderness, and Poisson's ratio on the natural frequencies of the rotor are investigated. Torabi et al 20 carried out the whirling analysis of axial-laded multistep Timoshenko rotor carrying concentrated masses, where forward and backward frequencies, and the corresponding mode shapes are obtained with the TMM. Afshari and Rahaghi 21 used the differential quadrature element method to study the transverse vibration of multistep Timoshenko rotors and the general conclusion was obtained: for a nonrotating rotor, the forward whirling speeds is equal to the backward ones, when the spin speed of the rotor increases, the forward ones increase and the backward ones decrease. AL-Shudeifat 22 found the new backward whirl zones of rotational speeds in intact and cracked rotor systems by using FEM. Afshari et al 23 studied the influence of lumped masses on the whirling speeds of a rotating Timoshenko shaft carrying multiple concentrated masses and found that: the influence of a concentrated mass on the forward and backward frequencies is nil if it is located at the "nodal points" of that natural mode shape and is highest when it is located at an "antinode" of that mode. By the way, the following References 24-28, may be useful for some readers. Huang et al 24  From the foregoing literature reviews, one sees that all the existing techniques for obtaining the forward or/and backward whirling speeds of the "shaft-disk systems" are either approximate or complicated except the simple analytical method presented in References 5,12. In theory, the solution given by Eshleman and Eubanks 5 is an exact one, however, it is only for the whirling speeds of a rotating shaft carrying "one" disk and the corresponding "whirling mode shapes" are not presented. Thus, in Reference 12, the above technique are extended and modified to obtain the lowest five (or higher) forward and backward whirling speeds and the corresponding mode shapes for a shaft carrying any number of disks with various (boundary) supporting conditions. Since the effects of SD, RI, and GM of the shaft are neglected in Reference 12, the purpose of this article aims at studying the title problem with SD, RI, and GM of the shaft considered. To this end, first of all, the mass moment of inertia for unit shaft length (at axial coordinate x), J (s) (x), and that for each rigid disk i, J D,i , are combined with their associated GMs to form the frequency-dependent equivalent mass moments of inertia, J (s) eq (x) and J eq, i , respectively. Next, the transverse displacement of each shaft cross-section is represented by a complex number, and then, for a rotating shaft-disk system, the equations of motion, the continuity equations for the deformations, the equilibrium equations for the forces (and moments), and the associated BCs are derived in terms of the complex numbers. In which, the effects of each rigid disk i is replaced by a lumped mass m d,i along with a frequency-dependent equivalent mass moment of inertia, J eq, i . After the foregoing manipulations, the whirling motion of a rotating three-dimensional (3D) shaft mounted by any number of rigid disks is similar to the transverse vibration of a stationary two-dimensional (2D) beam carrying any number of concentrated elements. Finally, the methods presented by the researchers [29][30][31][32][33][34] for the free vibration analyses of a stationary 2D beam carrying various concentrated elements can be used to determine the forward and backward whirling speeds and the associated whirling mode shapes of a loaded rotating 3D shaft. In addition to comparing with the existing literature, all results obtained from the presented method are also compared with those obtained from the FEM 35 and good agreements are achieved.
The numerical results obtained from the presented method are the exact solutions, they may be the benchmark for evaluating the accuracy of solutions of the other associated approximate methods (such as TMM, FEM, … ). In addition, the presented approach is available for the uniform or nonuniform shaft-disk systems. To reduce the deflections, most shafts are robust rather than slender, thus, the Timoshenko beam theory instead of the Euler beam theory should be used for the dynamic analyses of the title problem. Since the literature concerned is limited, the presented approach in this article should be valuable in the academic researches or the practical applications. Figure 1 shows a multistep shaft with its two ends supported by ball bearings. It is composed of n uniform shaft segments [denoted by (1), (2), … , (i−1), (i), (i+1), … , (n)], separated by n − 1 nodes [denoted by 1, 2, … , i−1, i, i+1, … , n−1] and carrying a rigid disk m d,i (with equivalent mass moment of inertia J eq, i ) at each node i, for i = 1 to n − 1. Figure 2 shows the coordinate systems for the rotating shaft with spin speed Ω about its longitudinal axis and whirling speedã bout the (horizontal) centerline of bearings (along x-axis). Where xyz,xỹz, and abc are the nonrotating (space-fixed), rotating (shaft-fixed), and cross-sectional (disk-fixed) coordinate systems, respectively. The relative position between xyz andxỹz coordinate systems is determined by the angle x =̃t, and the position of abc coordinate system is defined by x , u z / x, u y / x and the horizontal coordinate x, as one may see from Figure 2 of Reference 15. Besides, u y and u z are the vertical and horizontal transverse displacements of the cross-sectional centroid g of shaft (or center of gravity of disk) in y-and z-directions, while − u z / x and u y / x are the rotational angles of shaft (or disk) cross-section (located at x) about y-and z-axes, respectively.

F I G U R E 1 A multistep
bearing-support shaft composed of n uniform shaft segments [denoted by (1), (2), … , (i−1), (i), (i+1), … , (n)] separated by n − 1 nodes [denoted by 1, 2, … , i-1, i, i+1, … , n-1] and carrying a rigid disk m d,i (with equivalent mass moments of inertia J eq, i ) at each node i, for i = 1 to n − 1 F I G U R E 2 The coordinate systems for a shaft-disk system with spin speed Ω about a-axis and whirling speed about horizontal x-axis, with xyz,xỹz, and abc denoting the space-fixed, shaft-fixed, and disk-fixed (cross-sectional) coordinate systems, respectively From Equations (A21a,b) in the Appendix of this article, one obtains the equations of motion for the ith uniform spin Timoshenko shaft segment shown in Figure 1 to be In the above expressions, the primes ( ′ ) and the overhead dots (⋅) represent the differentiations w.r.t. the axial coordinate x and time t, respectively, while the symbols i , E i , G i , and A i are mass density, Young's modulus, shear modulus, and cross-sectional area of the ith shaft segment, respectively, I i is the diametric moment of inertia of area A i , u i (x,t) is the transverse displacement, i (x,t) is the rotational angle induced by the bending moment M i (x,t), i (x,t) is the shear strain induced by the shear force Q i (x,t), k ′ i is the shear coefficient and j = √ −1. For free vibrations, the functions, u i (x,t) and i (x,t), take the forms where U i (x) and Ψ i (x) are the amplitudes of u i (x,t) and i (x,t), respectively,̃is the whirling speed (rad/s) of the entire Timoshenko shaft and j = √ −1. Furthermore, in Equations (4a,b), the upper sign (+) and lower sign (−) are for forward and backward whirlings, respectively.
Substituting Equations (4a,b) into Equations (1a,b) and performing some manipulations, respectively, one obtains Equations (5a,b) take the same forms as the equations of motion for the stationary 2D bare Timoshenko beam, 32 and in Equation (6a), the symbol J (s) eq,i denotes the equivalent mass moment of inertia for unit length of the ith spin "shaft segment," it is similar to that for the ith rotating "disk," J eq, i , given later by Equation (21a).
For a circular shaft segment i with diameter d s,i , the moments of inertia for its cross-sectional area are given by and the associated mass moments of inertia per unit length are where J (s) i is the diametric mass moment of inertia per unit length and J (s) p,i is the polar one. Let i =̃2 then, the solutions of Equations (5a,b) take the forms The four integration constants (A i , B i , C i , and D i ) for the translational displacement function U i (x) given by Equation (8a) and those (A where For convenience, the displacement functions for translations and rotations given by Equations (8a,b) are rewritten below Furthermore, substituting the expressions for u i (x,t) and i (x,t) given by Equations (4a,b) along with that for i (x,t) given later by Equation (14e) into Equation (2b), and then inserting the expressions for U i (x) and Ψ i (x) given by Equations (12a,b) into the resulting expression, one obtains the displacement function for the SDs to be and, for free vibrations, one has introducing Equations (4b) and (14c,d,e) into Equations (14a,b) leads to , and i (x,t), respectively. For the first shaft segment with i = 1 and the final shaft segment with i = n (cf. Figure 1), the determination of the integration constants appearing in Equations (12a,b,c) (A i , B i , C i , and D i ) requires the following BCs: It is evident that the BCs of a shaft-disk system with various supporting conditions can be obtained from Equations (16a,b), (17a,b), and/or (18a,b), respectively.

CONTINUITY OF DEFORMATIONS AND EQUILIBRIUM OF FORCES (AND MOMENTS) AT EACH INTERMEDIATE NODE i OF A STATIONARY SHAFT
The continuity of displacements and slopes for the two shaft segments, (i) and (i + 1), joined at node i (see Figure 1) requires that Since a rotating 3D shaft-disk system can be replaced by an equivalent stationary 2D one, if each rigid disk i attached to node i of the spin shaft is replaced by an equivalent disk with lumped mass m d,i and frequency-dependent equivalent mass moment of inertia J eq, i, 12 the equations for the equilibrium of shear forces and bending moments of the two shaft segments, (i) and (i + 1), joined at disk i of the equivalent stationary 2D shaft-disk system are given by In Equation (21a), the upper sign (−) is for the forward whirling and the lower sign (+) is for the backward whirling. 12 For a circular thin disk with mass density d , diameter d d , and thickness h, its mass m d , diametrical mass moment of inertia J D and polar mass moment of inertia J p are given by Wu et al 12 Similarly, substituting U(x), Ψ(x), and Γ(x) defined by Equations (12a-c) into Equations (20a,b) leads to

BOUNDARY EQUATIONS AT TWO ENDS OF THE ENTIRE SHAFT
For the shaft-disk system shown in Figure 1, if it is supported by the rigid ball bearings, then the transverse displacements and bending moments at nodes 0 and n are equal to zero, thus, from Equations (18a,b) one has Introducing Equations (12a,b) into the above equations, one has A n cosh n L + B n sinh n L + C n cos n L + D n sin n L = 0 (30a) a n n A n cosh n L + a n n B n sinh n L − b n n C n cos n L − b n n D n sin n L = 0 (30b) Since the supporting condition of the current rotating 3D shaft is similar to that of the pinned-pinned (P-P) stationary 2D beam, it is called the P-P shaft in this article, for convenience.

DETERMINATIONS OF WHIRLING SPEEDS AND MODE SHAPES
For the shaft-disk system shown in Figure 1, the total number of equations, n, is equal to that of the integration constants of all shaft segments (A i , B i , C i , and D i , with i = 1 to n) given by n = 4n, and the simultaneous equations for all the integration constants take the form where {B} n×1 is a column vector composed of the n integration constants of all the n shaft segments, that is, and [H] n×n is a n × n square matrix.
Since the two ends of the shafting system shown in Figure 1 are supported by ball bearings (with the P-P BCs), the nonzero coefficients of [H] n×n may be obtained from Equations (29a,b), (23a,b), (24a,b), and (30a,b), respectively, as follows: (i) For the left pinned end at node 0 From Equations (29a,b) one obtains H 4n,4n−3 = a n n cosh n L, H 4n,4n−2 = a n n sinh n L, H 4n,4n−1 = −b n n cos n L, H 4n,4n = −b n n sin n L (40a-d) Now, the whirling speeds and the associated mode shapes of the shaft-disk system can be determined as follows: nontrivial solution for Equation (31) requires that The above expression is an eigenvalue equation, from which one may determine the whirling speeds of the shaft-disk system,̃r (r = 1, 2, 3, … ), by using the conventional half-interval method 36 or the modified one, 31 and corresponding to each whirling speed,̃r, one may obtain the associated integration constants, A i , B i , C i , and D i (i = 1 to n), from Equation (31). The substitution of the latter constants into Equation (12a) will define the corresponding rth whirling mode shape of the entire shaft,Ũ r (x) = ∑ n i=1Ũ r,i (x). For convenience, the equivalent mass moment of inertia for unit length of the ith spin "shaft segment" given by Equation . Similarly, the equivalent mass moment of inertia for each ith "disk" given by Equation (21a) is rewritten as . From the above two expressions one sees that the symbol C fb = ∓ 1.0, with upper sign "−" for the forward whirling motions and lower sign "+" for the backward ones. Thus, when one uses the half-interval method to find the lowest five "forward" whirling speeds̃F r (r = 1 − 5), one must set C fb = − 1.0 and then determine the lowest five "positive" roots of the equation | H(̃) |= 0 given by Equation (41). Similarly, when one determines the lowest five "backward" whirling speeds̃B r (r = 1 − 5), one must set C fb = + 1.0 and then determine the lowest five "positive" roots of the equation | H(̃) |= 0. It is noted that, in theory, one has "positive" natural frequencies and does not have the "negative" ones. In Equations (4a,b), the ± signs in e ±j̃t denote the (conventional) clockwise and counterclockwise directions of shaft whirling and have nothing to do with the magnitudes of̃. Furthermore, in the FEM, the same technique is used to find the values of F r (r = 1 − 5) and̃B r (r = 1 − 5). First, one sets C fb = − 1.0 and determines the values of̃F r (r = 1 − 5), next, one sets C fb = + 1.0 and determines those of̃B r (r = 1 − 5), in the above processes, all "negative" eigenvalues must be excluded.

NUMERICAL RESULTS AND DISCUSSIONS
In this section, the reliability of the presented theory and the developed computer programs for this article is confirmed first and then the characteristics regarding the stationary and the rotating uniform and nonuniform shafts carrying arbitrary disks are studied.

Comparisons with results of existing literature and FEM
If a shaft without any attachments is called the bare shaft and that carrying any rigid disks is called the loaded shaft, then, the former is the special case of the latter, because one may obtained the characteristics of the former from those of the latter by letting the mass density of each rigid disk to approach zero. Based on the last theory, first of all, the whirling speeds and mode shapes of a uniform P-P spin bare shaft are determined by using the presented analytical method and then compared with the existing literature to confirm the reliability of the presented theory and the developed computer programs. The data regarding the bare shaft are as follows 9 : The diameter d s of the shaft is not defined among the above expressions, because it is dependent on the Rayleigh beam coefficient C R = r g /L with r g denoting the radius of gyration of the shaft cross-section. The other parameter given by Zu and Han 9 is the nondimensional rotating speed defined by Ω = Ω∕ 10 , with Ω denoting the spin speed of the shaft and 10 denoting the first natural frequency of the corresponding stationary Timoshenko shaft.
Based on the data given by Equations (42a-f), with C R = 0.15 and L = 1.0, one obtains the radius of gyration of the shaft to be For a circular shaft with diameter d s , its radius of gyration r g is determined by If E = 207 Gpa = 2.07 × 10 11 N/m 2 and G = 77.6 Gpa = 0.776 × 10 11 N/m 2 , then the Poisson's ratio corresponding to For the case of spin speed Ω = 0, from column 2 of Table 1 one sees that the first natural frequency of the current stationary Timoshenko bare shaft is given by 10 Table 1. A comparison between the results of Table 1 and those of Reference 9 is shown in Figure 3, where the dotted curves ( … … ) are for the forward speed ratios  ãF r and̃B r refer to the rth "forward" and "backward" whirling speeds, respectively. b 10 = 2342.7942 rad/s is the fundamental natural frequency of the stationary P-P Timoshenko bare shaft.

F I G U R E 3 Influence of spin speed ratios
Ω/ 10 on the whirling speed ratios̃r∕ 10 of the P-P spin Timoshenko shaft, with dotted lines ( … … ) for first, second, or third "forward" speed ratios̃F r ∕ 10 and solid lines (--) for the "backward" ones̃B r ∕ 10 , all curves being obtained from Reference 9. Besides, the symbols, •, +, and ▲, for the first, second, and third "forward" whirling speed ratios, and those, ○, ×, and △, for the first, second, and third "backward" ones being obtained from this article B r =̃B r ∕ 10 (r = 1, 2, 3), and all the above curves are obtained from Figure 2 of Reference 9. In addition, the symbols • and ○ denote the "first" forward and backward whirling speed ratios, those + and × denote the "second" ones, and those ▲ and △ denote the "third" ones, respectively, and all the last symbols are obtained from the results of this article (as shown in Table 1). From Figure 3 one sees that the results of this article are in good agreement with the corresponding ones of Reference 9. F I G U R E 4 The lowest three whirling mode shapes of the spin Timoshenko shaft,Ũ r (x∕L) (r = 1, 2, 3), corresponding to the nondimensional whirling speeds̃r∕ 10 (r = 1, 2, 3) shown in Figure 3 and Table 1  For the case of spin speed ratio Ω = Ω∕ 10 = 2.5, the lowest three forward and backward unit-amplitude mode shapes, U F r (x∕L) andŨ B r (x∕L) (r = 1, 2, 3), of the bare shaft corresponding to the nondimensional whirling speeds̃r∕ 10 (r = 1, 2, 3) shown in Figure 3 and Table 1, are plotted in Figure 4. In which, the legends are similar to those for Figure 3, that is, the dotted and solid curves are obtained from Figure 4 of Reference 9 and the symbols are obtained from the results of this article. It is seen that each pair of forward and backward whirling mode shapes overlap each other, and this result agrees with Reference 9.
Since the main difference between the whirling motions of a Timoshenko shaft-disk system (with SD, RI, and GM of the shaft considered) and those of an Euler shaft-disk system (with SD, RI, and GM of the shaft neglected) is due to the spin shaft, thus, the good agreement between the lowest three whirling speeds,̃F r and̃B r (r = 1, 2, 3), as well as the associated forward and backward mode shapes,Ũ F r (x∕L) andŨ B r (x∕L) (r = 1, 2, 3), of the spin bare Timoshenko shaft obtained from this article and those obtained from Reference 9, as shown in Table 1 and Figures 3 and 4, should be one of the good evidences for the reliability of the results presented in this article.
In the last example, only the characteristics of a bare shaft are investigated, and those of a loaded shaft will be studied in this example. Figure 5 shows the mathematical model of the shaft-disk system considered. The given data are as follows 12 :  In the above Equations (46a-h), the subscripts "s" and "d" refer to "shaft" and "disk," respectively. In addition, the shear modulus of shaft material is determined by (cf. Equation (44) The data given by Equations (46a-h) are the same as those for Table 1 of Reference 12, except that the shaft studied in the last reference is an Euler beam, so that the effects of SD, RI, and GM of the spin shaft are neglected. 12 Based on the Timoshenko shaft theory presented in this article and the data given by Equations (46) and (47), the lowest five natural frequencies ( 1 − 5 ) of the stationary shaft-disk system (with = 0) are listed in Table 2A and the lowest five whirling speeds (̃1 −̃5) with = 1.0 are listed in Table 2B. From Table 2A one sees that the lowest five natural frequencies obtained from the presented (exact) method are very close to those obtained from FEM, furthermore, because of the effects of SD and RI, the values of 1 − 5 are smaller than the corresponding ones obtained from the Euler shaft-disk system 12 as shown in the final row of Table 2A. For the spin Timoshenko shaft-disk system with = 1.0, from Table 2B one sees that, either forward whiling speeds (̃F 1 −̃F 5 ) or backward ones (̃B 1 −̃B 5 ), the results obtained from the presented method are very close to the corresponding ones obtained from FEM. In addition, the whirling speeds of the spin Timoshenko shaft-disk system are also smaller than the corresponding ones of the spin Euler shaft-disk system 12 shown in the final two rows of Table 2B due to the effects of SD and RI. From Table 2A,B, one finds the relationships̃F 1 =̃B 1 = 1 and F 3 =̃B 3 = 3 for the spin Euler shaft-disk system as mentioned by Wu et al, 12 however, they are not true for the spin Timoshenko shaft-disk system studied in this article except that the GM of the spin bare Timoshenko shaft is neglected. This is a reasonable result, because the slopes of the first and third natural mode shapes at x = x 1 = L/2 (where the single rigid disk located) are no longer equal to zero due to the effect of GM of the spin Timoshenko shaft. Based on the foregoing results, it is believed that the formulations and the developed computer programs associated with the presented (exact) method and the FEM should be reliable.  Table 2B, one sees that the whirling speeds obtained from the presented (exact) method,̃F r,Exact (r = 1 − 5), listed in the third row for the case with GM ≠ 0 are greater than the corresponding ones listed in the seventh row for the case with GM = 0, this result indicates that the GM of the spin Timoshenko shaft has the effect of raising its forward whirling speeds̃F r,Exact (r = 1 − 5). On the contrary, from the values of̃B r,Exact (r = 1 − 5) listed in the fifth row of Table 2B for the case with GM ≠ 0 and the corresponding ones listed in the 9th row for the case with GM = 0, one sees that the GM has the effect of reducing the backward whirling speeds̃B r,Exact (r = 1 − 5).
Corresponding to the lowest five natural frequencies listed in Table 2A with = 0 and the lowest five whirling speeds listed in Table 2B with = 1, the associated natural mode shapes and whirling ones of the P-P Timoshenko shaft-disk system are shown in Figures 6 and 7, respectively. From Figure 6 one sees that the lowest five unit-amplitude mode shapes U 1 ( ) − U 5 ( ) (with = x/L) obtained from the presented method and denoted by the solid lines (--) are very close to the corresponding ones obtained from FEM and denoted by the dotted lines ( … … ). In which, the first, second, third, fourth, and fifth natural mode shapes are denoted by the symbols •(or ○), +(or ×), ▲(or △), ◼(or □), and ★ (or ⋆), respectively. Similarly, either the lowest five forward whirling mode shapesŨ F 1 ( ) −Ũ F 5 ( ) shown in Figure 7A or the lowest five backward onesŨ B 1 ( ) −Ũ B 5 ( ) shown in Figure 7B, those obtained from the presented method (denoted by the solid lines) are very close to the corresponding ones obtained from FEM (denoted by the dotted lines). The legends for Figure 7 are the same as those for Figure 6.

Whirling motions for a uniform spin shaft carrying three identical disks
In this section, the natural frequencies and associated mode shapes of a stationary P-P uniform shaft (with = 0) carrying three identical rigid disks (each with thickness h = 0.012 m) as shown in Figure 8 are studied first and shown in Table 3A and Figure 9, then, the whirling speeds and associated mode shapes of the spin shaft-disk system (with = 1) are investigated and shown in Table 3B and Figure 10. Next, the influence of spin speeds (Ω) of the shaft on the rth whirling speeds (̃r) of the entire shaft-disk system are determined and shown in Figure 11. All data for the present example are the same as those for the last example given by Equations (46) and (47), except that the diameter of the shaft is larger, that is, d s = 0.03 m (instead of d s = 0.02 m). Unless particularly mentioned, hereafter, all shafts refer to the Timoshenko shafts with effects of SD, RI, and GM considered.  Table 3A (with = 0) one sees that the lowest five natural frequencies of the stationary shaft-disk system (Figure 8) obtained from the presented method are very close to those obtained from the FEM, and this is the reason why the associated mode shapes obtained from the presented method (the solid lines) and those obtained from the FEM (the dashed lines) overlap each other in Figure 9. Besides, from Table 3B, one sees that, for the lowest five forward whirling speeds,̃F r,Exact (r = 1 − 5), obtained from the presented (exact) method and listed in the third and seventh rows, those  listed in the third row with the GM of shaft considered are greater than the corresponding ones listed in the seventh row with the GM of shaft neglected, this indicates that the GM of the shaft has the effect of "raising" the forward whirling speeds. However, the last effect of GM is reverse for the backward ones̃B r,Exact (r = 1 − 5), because the lowest five backward whirling speeds listed in the fifth row with the GM of the shaft considered are "smaller" than the corresponding ones listed in the ninth row with the GM of the shaft neglected. It is also noted that, in Table 3A,B, the results obtained from the presented method are very close to the corresponding ones obtained from FEM. Corresponding to the lowest five whirling speeds listed in the third and fifth rows of Table 3B, the associated "forward" whirling mode shapesŨ F r ( ) with = x/L and "backward" onesŨ B r ( ) for the spin loaded shaft are plotted in Figure 10. The meanings for the symbols are the same as those appearing in Figures 6, 7, and 9, that is, the first, second, third, fourth, and fifth whirling mode shapes are denoted by the symbols •(or ○), +(or ×), ▲(or △), ◼(or □), and ★ (or ⋆), respectively. It is seen that, except for the first mode shapes (Ũ F 1 ( ) andŨ B 1 ( )), for the other higher ones, the mode shapesŨ F r ( ) of "forward" whirling are different from the corresponding onesŨ B r ( ) of "backward" whirling, for r = 2 − 5. The last result agrees with that obtained from the spin Euler shaft-disk system. 12 Furthermore, by means of the presented method, the influence of spin speeds Ω on the lowest five "forward" whirling speeds̃F r (r = 1 − 5) represented by the solid lines (with symbols •, +, ▲, ◼, and ★) and the "backward" ones̃B r (r = 1 − 5) represented by the dashed lines (with symbols ○, ×, △, □, and ⋆) is shown in Figure 11. It is seen that the influence of spin speeds Ω on the whirling speeds of the loaded F I G U R E 9 The lowest five natural mode shapes of the stationary uniform P-P shaft ( = 0) carrying three identical rigid disks (Figure 8) obtained from presented method (solid lines with symbols •, +, ▲, ◼, and ★) and finite element method (dashed lines with symbols ○, × , △, □, and ⋆)

F I G U R E 10 The lowest five
forward whirling mode shapes (solid lines) and backward ones (dashed lines) for the spin uniform P-P shaft carrying three rigid disks (Figure 8 Figure 11 is similar to that of the spin bare shaft shown in Figure 3. This is a reasonable result, because both the loaded shaft and the bare shaft are symmetrical with respect to their central cross-sections (located at x = L/2). It is seen that, in Figure 11, the spin speed Ω corresponding to the intersection of the rth-mode curves (for̃F r and B r ) and the (vertical) ordinate (axis) is zero (ie, Ω = 0), thus, the associated whirling frequencies are equal to the natural frequency, that is,̃F r =̃B r = r at Ω = 0. In addition, for each pair of rth-mode curves, the curve for forward speeds F I G U R E 11 The influence of spin speeds Ω on the lowest five forward and backward whirling speeds (̃F r and̃B r , r = 1 − 5) of the spin uniform P-P symmetrical shaft-disk system (Figure 8) obtained from the presented method F r and that for backward speeds̃B r are approximately symmetrical with respect to the associated horizontal dotted line passing through the point on the ordinate (axis) with̃F r =̃B r = r at Ω = 0. This indicates that the influence of spin speeds on the forward motions is the same as that on the backward motions for a "symmetrical" shaft-disk system like that shown in Figure 8. It is noted that the "unit" for̃F r and̃B r (r = 1 − 5) in Figures 11 and 15 is "Hz" and 1 Hz =2 rad/s. The plot like Figure 11 to show the influence of rotating speeds Ω (unit Hz) on the rth whirling speeds̃X r (unit Hz), with X = F for forward whirling and X = B for backward whirling, is called the Campbell diagram. 19,20 In which, any line with slope equal to unit (ie, Δ̃X r ∕ = 1.0) is called the synchronous whirling line. 19,20 It is noted that the intersection of the last line with the Campbell diagram determines the critical speeds which should be avoided. For the Campbell diagram shown in Figure 11, the coordinates at the origin are O(0, 0) and those for the opposite point are P (300, 300), thus, the red thick dashed line connecting points O and P in Figure 11 represents a "reference" synchronous whiling line with slope Δ̃X r ∕ = (300 − 0)∕(300 − 0) = 1.0. It is evident that any line parallel to the last red thick dashed line in Figure 11 is a synchronous whirling line.

Whirling motions for a nonuniform spin shaft carrying three identical disks
To show the effectiveness of the presented method for solving the whirling problems of a nonuniform shaft carrying arbitrary disks, a two-step shaft carrying three identical disks shown in Figure 12 is studied in this section. All data are the same as those of the last example ( Figure 8) except that the diameter of shaft segment from x = 1 3 L ∼ 2 3 L is replaced  by the one with larger diameter d s2 = 0.04 m (instead of d s1 = 0.03 m). In addition, the locations for the first and final disks are also slightly changed. For the bare shaft, from Table 4A and Figure 13 one sees that the lowest five natural frequencies r (r = 1 − 5) and the corresponding mode shapes U r ( ) (r = 1 − 5) obtained from the presented method are very close to those obtained from FEM. Furthermore, the lowest five natural frequencies listed in Table 4A for the nonuniform shaft-disk system shown in Figure 12 are higher than the corresponding ones listed in Table 3A for the uniform shaft-disk system shown in Figure 8, because the stiffness of the former is greater than the latter. However, the mode displacements of the third and fifth mode shapes near the center (at x = L/2) of the nonuniform shaft shown in Figure 13 are smaller than those of the uniform shaft shown in Figure 9, this is because a bigger shaft segment locates at x = 1 3 L ∼ 2 3 L in the nonuniform shaft. For the loaded shaft, from the lowest five forward whirling speeds̃F r (r = 1 − 5) and backward ones̃B r (r = 1 − 5) listed in Table 4B, one sees that the previous statements about the whirling speeds of the uniform shaft-disk system shown in Table 3B are also available for the current nonuniform shaft-disk system. For example, for the forward whirling, the lowest five speeds̃F r (r = 1 − 5) with GM of the shaft "considered" and obtained from either presented method (or FEM) and listed in the third (or fourth) row of Table 4B are greater than the corresponding ones with GM of the shaft "neglected" and listed in the seventh (or eighth) row of Table 4B. On the other hand, for the backward whirling, the lowest five speeds F r (r = 1 − 5) with GM of the shaft "considered" are smaller than the corresponding ones with the GM of the shaft is F I G U R E 13 The lowest five natural mode shapes of the stationary nonuniform (two-step) P-P shaft carrying three identical rigid disks (see Figure 12) obtained from the presented method. The legends are the same as those for Figure 9 "neglected," as one may see from the fifth (or sixth) and the 9th (or 10th) row of Table 4B. It is noted that the average CPU time required by the presented method is approximately equal to one half of that required by FEM (by using an IBM PC Pentium III) as one may see from Table 4A,B. In addition, from the whirling mode shapes of the nonuniform shaft-disk system plotted in Figure 14, one sees that, the forward mode shapesŨ F r ( ) (r = 1 − 5) are different from the backward ones U B r ( ) (r = 1 − 5) except for the first mode shapesŨ F 1 ( ) andŨ B 1 ( ). Finally, from the curves for the lowest five forward and backward whirling speeds (̃F r and̃B r , r = 1 − 5) vs the spin speeds Ω shown in Figure 15 obtained from presented method, one sees that the influence of spin speeds on the forward whirling speeds is the same as that on the backward ones. It is also noted that the curve for the fourth forward whirling speeds̃F 4 intersects the curve for the fifth backward whirling speeds̃B 5 at Ω ≈ 120 Hz. This indicates that the fourth forward-mode response and the fifth backward-mode response of the shaft-disk system may be coupled each other if Ω ≥ 120 Hz. It is similar to Figure 11 that, in Figure 15, the red thick dashed line connecting the origin O(0, 0) and the opposite point P(300, 300) represents a "reference" synchronous whirling line for the Campbell diagram and each line parallel to the red thick dashed line is a synchronous whirling line.

Influence of disk positions on whirling speeds for a uniform shaft carrying one disk
Since the whirling speeds (̃F r and̃B r with r = 1, 2, 3, … ) of a spin loaded Timoshenko shaft carrying one disk is dependent on the magnitude of the GM developed by the attached disk, and the latter is dependent on the "slope" U ′ r (x) of the mode shape where the disk is located, 12 the "key positions" for the disk should be at the points on a mode shape (of the shaft-disk system) with "slope" to be minimum or maximum. Based on the above statements and the lowest five natural mode shapes of the stationary "bare" P-P Timoshenko shaft shown in Figure 16, we select the next positions for the disk: = x/L = 0.25 and 0.5, and the reasons are as follows. From Figure 16 one sees that at = x/L = 0.25, the slope of second mode shape approaches zero (ie, U ′ 2 (0.25L) ≈ 0) and that of third mode shape approaches 1.0 (ie, | U ′ 3 (0.25L) |≈ 1.0). Furthermore, at = x/L = 0.5, the slope of first or third mode shape approaches zero (ie,U ′ r (0.5L) ≈ 0, r = 1, 3) and that of second one approaches 1.0 (ie, | U ′ 2 (0.5L) |≈ 1.0). The foregoing mode shapes are for the stationary "bare" P-P Timoshenko shaft, and the mode shapes for the stationary "loaded" one are shown in Figures 17 and 6 for the cases of rigid disk located at F I G U R E 14 The lowest five forward whirling mode shapes and backward ones for the spin nonuniform P-P shaft carrying three rigid disks (see Figure 12) obtained from the presented method with = 1.0. The legends are the same as those for Figure 10 Figure 16. However, the positions with slope to be near 0 or 1.0 are not significantly changed, particularly for those appearing in Figure 6. It is noted that, in this subsection, all input data are the same as those given by Equations (46) and (47).
By using the presented method, the influence of disk positions on the lowest five natural frequencies, 1 − 5 , for the uniform "stationary" loaded P-P shaft (with = 0) carrying one rigid disk at = x/L = 0.25 and 0.5 (cf. Figure 5) are listed in Table 5A, and for the case of = 1.0, the influence of disk positions on the lowest five whirling speeds,̃X 1 −̃X 5 (X = F,B) is shown in Table 5B. For the case of = x/L = 0.25, from Figure 17 and Table 5B, with respect to the five TA B L E 5 Influence of disk positions (ξ = x/L) on the lowest five (A) natural frequencies (with λ = 0), ω 1 − ω 5 , and (B) whirling speeds (with λ = 1.0),̃1 −̃5, for the uniform loaded P-P shaft carrying one rigid disk at ξ = x/L (cf. Figure 5)  natural frequencies listed in Table 5A, one sees that the forward whirling speeds increase gradually and the trend for the backward ones is reverse. On the other hand, for the case of = x/L = 0.5, from Figure 6 and Table 5B, with respect to the five natural frequencies listed in Table 5A, one sees that̃F 1 ≈̃B 1 ≈ 1 and̃F 3 ≈̃B 3 ≈ 3 , and the difference of̃F 2 −̃B 2 increases significantly. To confirm the foregoing results obtained from Figures 17 and 6 and Table 5A,B, the influence of spin speeds Ω on the lowest four forward and backward whirling speeds (̃F r and̃B r with r = 1 − 4) of the spin uniform P-P shaft-disk system (cf. Figure 5) obtained from the presented method are shown in Figure 18A for the case of = x/L = 0.25 and Figure 18B for = x/L = 0.5, respectively. It is easy to see that all conclusions drawn from Figures 17 and 6 and Table 5A,B are satisfied by the curves shown in Figure 18A,B.

CONCLUSIONS
1. For a three-dimensional "rotating" shaft-disk system (or simply called 3D shaft-disk system), if the GM of shaft is combined with its diametric mass moment of inertia (J (s) ) to form an equivalent mass moment of inertia of the shaft per unit length (J (s) eq ), and the GM of each disk i is also combined with mass moment of inertia of each disk i (J D,i ) to form an equivalent mass moment of inertia of each disk i (J eq, i ), then the 3D shaft-disk system can be replaced by a 2D "stationary" beam carrying the same number of disk(s) so that the whirling speeds and the associated whirling mode shapes of the 3D shaft-disk system can be easily obtained. 2. The formulation for the whirling motion of a spin "Timoshenko shaft-disk system" presented in this article is also available for the dynamic analysis of a bare uniform or nonuniform spin "Timoshenko shaft," the only requirement is to set the mass density of each "disk" mounted on the Timoshenko shaft to approach zero. 3. Due to the effects of SD and RI of the shaft, the lowest five natural frequencies r (r = 1 − 5) of the stationary Timoshenko shaft-disk system are smaller than the corresponding ones of the associated stationary Euler shaft-disk system. For the same reason, the lowest five forward (or backward) whirling speeds̃X r (X = F,B and r = 1 − 5) of the spin Timoshenko shaft-disk system are also smaller than the corresponding ones of the associated spin Euler shaft-disk system no matter whether or not the GM of the shaft is considered. 4. For a shaft-disk system with spin speed Ω, the GM of either its shaft or its disk(s) can raise the forward whirling speeds, F r (r = 1 − 5), and reduce the backward ones,̃B r (r = 1 − 5), this is the reason why the curve of̃F r vs Ω intersects the curve of̃B r vs Ω at the point on the ordinate (axis) with Ω = 0 for each specified rth mode, and each pair of curves is approximately symmetrical with respect to the horizontal line passing through the point with̃F r =̃B r = r (and Ω = 0) for a "symmetrical" Timoshenko or Euler shaft-disk system. 5. In theory, the numerical results obtained from an analytical method are the exact solutions, thus, the solutions of the analytical method presented in this article may be the benchmark for evaluating the accuracy of the other associated approximate methods (such as FEM, TMM, … ). Furthermore, the presented method is available for the determination of accurate whirling speeds and the associated mode shapes of a uniform or nonuniform rotating shaft carrying arbitrary disks in various BCs. 6. For a shaft-disk system with spin speed Ω, the rth forward whirling speed̃F r and the rth backward onẽB r , if thẽF r vs Ω curve intersects thẽB r+1 vs Ω curve at Ω = Ω*, then the rth forward-mode response and the (r+1)th backward-mode response may couple each other when Ω ≥ Ω* 7. The influence of position for a rigid disk in a shaft-disk system is dependent on the "slope" of the rth (shaft) mode shape at the point where the rigid disk is located, the larger the "slope" the larger the influence of the spin shaft on the 30. Lin HY. On the natural frequencies and mode shapes of a multi-span and multi-step beam carrying a number of concentrated elements. How to cite this article: Wu J-S, Hsu T-F. An efficient approach for whirling speeds and mode shapes of uniform and nonuniform Timoshenko shafts mounted by arbitrary rigid disks. Engineering Reports. 2020;2:e12183. https://doi.org/10.1002/eng2.12183 Figure A1 shows the coordinate systems for a rotating bare shaft with spin speed Ω about its longitudinal axis and whirling speed̃about the horizontal (x) axis along the bearing centerline. Where oxyz is a nonrotating (space-fixed) coordinate system,õxỹz is a rotating (shaft-fixed) coordinate system and gabc is a rotating coordinate system located at the shaft cross-section. Besides, u y and u z are the transverse displacements for the centroid g of the shaft cross-section (located at x) in y-and z-directions, respectively. It is noted that all coordinate systems shown in Figure A1 are based on the right-hand screw rule, thus, the positive rotational angle z on the xy-plane is about +z-axis, but the positive rotational angle y on the xz-plane is about −y-axis.

APPENDIX A.. EQUATIONS OF MOTION FOR A UNIFORM SPIN BARE TIMOSHENKO SHAFT
With respect to the stationary oxyz coordinate system, if the deformations of the shaft are small, then referring to Figure A2A,B, one may obtain the kinetic energy T, the strain energy V and the nonconservative virtual work W nc of the entire shaft given by Han and Zu. 10 W nc = f y u y + f z u z + m y y + m z z (A3)

F I G U R E A1 A rotating bare shaft
with spin speed Ω about its longitudinal a-axis and whirling speed̃about the horizontal x-axis, where oxyz is the nonrotating (space-fixed) coordinate system andõxỹz is the rotating (shaft-fixed) coordinate system. Besides, gabc is a rotating (section-fixed) coordinate system located at the shaft cross-section (with horizontal coordinate x) For free vibrations, f y = f z = 0 and m y = m z = 0, thus, Equations (A5a,b) and (A6a,b) become k ′ GA (u ′′ y − ′ z ) − Aü y = 0 ( A 9 a ) k ′ GA (u ′′ z − ′ y ) − Aü z = 0 (A9b) For a circular shaft with diameter d, the values of I y , I z , and I x are given by  Differentiating Equation (A16) w.r.t. x, substituting the expression of u ′ given by Equation (A20b) into the resulting expression and then using the expression of given by Equation (A20c), one obtains For convenience, Equation (A17) is rewritten below