Modelling of rotors with axisymmetric solid harmonic elements

Journal ofSound and Vibration (1989) 131(3), 431-443

MODELLING OF ROTORS WITH AXISYMMETRIC HARMONIC ELEMENTS R. W.

STEPHENSON,

SOLID

K. E. ROUCH AND R. ARORA

Department of Mechanical Engineering, University of Kentucky, Lexington, Kentucky 40506, U.S.A. (Received 25 July 1988, and in revised form 22 November 1988)

In rotating machinery analysis, the rotating shaft is typically modeled by a series of line or beam elements. These beam elements are formulated from classical beam theory, with a basic assumption that “plane sections remain plane” during bending. However, when the rotor has abrupt changes in the diameter, this assumption is not valid; local distortion occurs near the section change, resulting in an increase in the overall bending flexibility of the shaft. In this paper the bending stiffness and natural frequencies of several shaft geometries with stepped diameter changes are compared using several modelling assumptions. The comparison includes beam approximations with both transfer matrix and finite element approaches, and available experimental data. In addition, the conceptual approach of using an axisymmetric solid model of the shaft is reviewed. Results are included for a commercial four-node element, and for a cubic element developed for this work. The effective use of matrix reduction in conjunction with an axisymmetric model is shown. The importance of non-planar distortion near section changes is reviewed.

1. INTRODUCTION

Accurate prediction of rotor system dynamic characteristics is essential for design of modern rotating machinery, and numerical approximations are commonly used for complex systems. The two usual numerical approaches are the transfer matrix method [I, 21 and the finite element method [3-61. Both have been used successfully. Typically, a shaft model (finite element or transfer matrix) consists of a series of line elements with length and diameter specified to represent each section of the shaft. A number of diff erent beam element formulations have been developed, most including the effect of shear and rotary inertia. Either lumped or consistent mass matrices can be used [7,8]. However, all of these are derived from classical beam theory, which specifies that plane sections remain plane during bending. This assumption, while valid for a uniform shaft, is not necessarily valid in the vicinity of a stepped diameter change. Distortion occurs in the cross-section such that the bending strain is no longer a linear function of the radius. This local distortion increases the overall bending flexibility of the shaft. Typical construction in modem turbomachinery includes a number of diameter changes to accommodate bearings, seals, and shrunk-on impellers, fans and disks. For accurate predictions of the shaft dynamic response it is necessary to include this local reduced stiffness effect. Otherwise, the calculated critical speeds and other dynamic characteristics may be in error. In this paper, the bending stiffness and bending natural frequencies of several shafts with stepped diameter changes are evaluated using (a) classical beam elements and (b) finite element results from models using axisymmetric solid (harmonic) elements which allow distortion of the cross-section. Comparisons are also made with experimental data. 431 0022-460x/89/120431+13

%03.00/O

@ 1989 Academic Press Limited

432

R. W. STEPHENSON

ET AL.

2. BACKGROUND The simplest beam formulations are based on the Bernoulli-Euler assumptions in which shear deflection and rotary inertia are not included. However, if a beam element becomes short relative to its diameter, or if higher vibration modes are involved, these effects become significant [9]. The Timoshenko beam theory includes shear and rotary inertia and forms the basis of most rotor dynamics analysis. The beam formulation can be in terms of either the transfer matrix or finite element method. By a transfer matrix, a relationship is established between a set of unknowns on the two ends of the beam section. The equations for the entire system are obtained by multiplying the successive matrix expressions for connected elements. Transfer matrix methods have been extensively used since the 194Os, with the Holzer-Myklestad and Prohl work. The finite element method has developed in parallel with the expanded computational power during the past two decades. In finite element, a stiffness (and mass and damping) matrix relate a displacement (and time derivatives) to forces and moments for degrees of freedom on each element. The system matrix representation is obtained by adding element matrices at appropriate entries dictated by the connectivity. Again, the individual element matrices reflect the assumptions involved in their formulation. For the finite element method, there are a number of beam formulations, primarily based on the Timoshenko assumptions. A review of these has been given by Thomas, Wilson and Wilson [lo]. More recently, tapered beam elements have been developed by Rouch and Kao [4] and Greenhill et al. [5] to model a linearly varying diameter along the beam length. In the formulations, the shear deflection is typically treated by assuming a linear variation of strain versus position in a finite element by either retaining shear strain as a degree of freedom or incorporating its effect into the translational motion at the nodes. The effect of the shape cross-section on shear deflection is taken care of by a shear deflection factor. For example, the circular sections in this presentation were assigned a factor of 10/9. In most of the beam finite element formulations a consistent mass formulation is used. In the transfer matrix method, a lumped mass approach is often employed, such that both the mass and inertia properties are lumped at discrete points along the beam, which are then connected by massless elastic members. A lumped mass approach is more convenient for computational processes. However, one advantage of a consistent mass approach lies in its convergence characteristics. In natural frequency analyses of a conforming element, the calculated frequencies will always be above the true frequencies, and will converge from above to the true frequencies as the discretization of the rotor is increased. This convergence from above characteristic does not exist with a lumped mass approach (transfer matrix or finite element). For any of the beam element formulations, a fundamental assumption is made that “plane sections remain plane during bending”; i.e., at a cross-section the axia1 strains vary in a linear manner as a function of the radial position. This is a valid assumption for a uniform circular shaft. In the vicinity of a stepped diameter change, however, a degree of local distortion occurs in the cross-section and the effect is that the local flexibility is increased. Studies on this effect for static deflection have resulted in curves such as those by Brown [ll], giving an effective diameter of the shaft based on the non-dimensional values of D/d and L/d, as shown in Figure 1. This effective shaft diameter, when used along the length of the large diameter section, will approximate the stiffness in a beam model of the shaft for rotor dynamic analysis. In the case of reference [ 111 the stepped increase is followed by an equal stepped decrease. Another possible geometry consists of a single step as shown in Figure 2. Here, the large diameter becomes

ROTOR

MODELLING

WITH

SOLID

Figure 1. Shaft correction factor with stepped diameter 0~036(L/d)[1~0-e~2’6Z5r’D~d’-‘1”L”()], where D, = Denec,,ve.

ELEMENTS

change,

433

from Brown [ 111, De/d = l.O+

Transition &on

Figure 2. Shaft with one local diameter change.

fully effective only after a certain distance from the step change. A reduced effective diameter could be used in this transition region, although no correction factors were found for this type of geometry. Rotor analysts often apply such corrections, or more approximate methods, such as the “45 degree rule”. An alternative technique presented in this paper uses axisymmetric solid harmonic elements, instead of beam elements, for modelling rotors with stepped diameter changes. This eliminates the need for correction factors since the cross-section distortion is modeled in detail during bending for distortion of the cross-section. In the axisymmetric elements discussed here, all displacements (and other quantities) are expressed as u(x, y, 6) = V(X,y) cos n6 or u(x, y, 0) = u(x, y) sin no, where 0 represents the angular co-ordinate in the circumferential direction so that a Fourier series representation is achieved [12]. In the case of a beam in lateral bending, the only mode of interest is that of n = 1, and the superposition of the series terms can be avoided. With this approach, the model is only required to represent the x-y plane, and the circumferential variation is cared for in the element itself. The 2D form of the element can be thought of as a complete ring in its mathematical representation. The three-dimensional problem is converted into a two-dimensional problem, with all of the geometry effects in the cross-section retained. For a node on the axis of revolution (zero radius), the axial degree of freedom needs to be constrained to zero for the lateral bending mode n = 1 so that the strains remain finite. For a free rotor, or one with equal stiffness in the orthogonal directions, only the symmetric radial, axial and circumferential degrees of freedom are required. In an asymmetric support situation, the antisymmetric degrees of freedom, as defined in reference [13], must also be included. Geradin and Kill [13] have used an axisymmetric element for rotor dynamics in conjunction with the classical component mode method of system reduction by splitting the rotor model into a number of substructures. An alternative to the component mode method for reducing the system dynamic degrees of freedom is the Guyan reduction technique [15,16]. In Guyan reduction, the system equations are in terms of a reduced number of degrees of freedom. The mass and stiffness of the condensed degrees of freedom, or slaves, are related to the retained degrees of freedom, or masters. Basic guidelines for choosing the master degrees of freedom are discussed in references 161 and [IS]. An

434

R. W. STEPHENSON

ET AL.

axisymmetric finite element for rotor dynamics, using Guyan reduction, has been developed [14]. The mass and stiffness matrices for the element in this program are summarized in the Appendix. While it is recognized that the solid element approach involves many more degrees of freedom than that of a beam approach, a reduced matrix approach can provide a much smaller system of equations without loss of accuracy in dynamic analysis [3,6]. The final reduced matrix size with the axisymmetric elements need not be significantly larger than a beam formulation.

3. RESULTS

3.1. STATIC STIFFNESS The static characteristics of a beam with varying cross-section will be examined initially, to prove a basis for discussion of dynamic characteristics. The finite element model is treated as a good approximation of the actual deflection in the comparison. The calculated stiffness values for a cantilever beam, shown in Figure 3, with a diameter change are given in Table 1 for D/d = 2.0 and L/d = 2-O (Figure l), as obtained by using beam and finite element models. For transverse loading, the error in the stiffness caused by neglecting the diameter change can be as significant as the error caused by neglecting shear deflection effects. For the load case with an applied moment at the end of the beam there is no shear deflection. The correction factor (Brown [ 111) significantly improves the accuracy of the bending stiffness for both load cases. The results from beam theory with correction factors differ slightly from the finite element results, apparently due to the approximations in the correction factor. For the finite element model with solid axisymmetric elements [ 151, the axial displacements through the cross-section at different axial locations can be examined. Figure 4 is a finite element model of the beam in Table 1, with axisymmetric solid elements. The

50.8++--50,8--+-50.8

Figure 3. Beam dimensions TABLE

(mm) used in

Table 1.

1

Static stiflness of cantilever beam with diameter change Stiffness

‘Transverse load (N/m)

Moment (Nm/rad)

Beam theory: (a) without shear, without correction factor (b) without shear, with correction factor (c) with shear, without correction factor (d) with shear, with correction factor Finite element,

axisymmetric

solid elements

4.73 E6 4.55 E6 4.66 E6 4.48 E6 4.38 E6

4.03 E4

3.81 E4 3.74 E4

ROTOR MODELLING

435

WITH SOLID ELEMENTS

Figure 4. Finite elementmodel of Figure 1 configuration.

end of the beam near point E is fixed and it has a concentrated bending moment applied at the other end which is free. Table 2 shows the axial displacements, as a function of radius at different axial locations (labeled in Figure 4) along the beam, resulting from the applied moment. The displacements are normalized with respect to the axial displacement at the outside diameter at each cross-section. At locations well removed from a diameter change, the displacements are linearly proportional to the radial location, confirming that plane sections remain plane in accordance with classical beam theory. However, in regions near a diameter change this is not true, indicating a warping of the cross-section. This additional flexibility lowers the stiffness locally. The effect on the total beam stiffness depends on the extent of this local flexibility relative to the overall bending flexibility of the beam. Figure 5 is a finite element model of the geometry in Figure 2 with D/d = 2.0. The axial displacements at different cross-sections are given in Table 3 for an applied moment TABLE 2

Axial displacements of beam in Figure 4 Axial location

Non-dimensional displacement

Beam theory displacement

0.333 0.667

0.350 0.700 1 .oo

0.333 0.667 1.00

5.1 4.9 0.0

0.333 0.667 1.67 2.00

0.188 0.372 0.539 0,677 0.838 1.00

0.167 0.333 0.500 0.667 0.833 1.00

12.5 11.7 7.8 1.5 0.6 0.0

C

0.333 0.667 1.00 1.33 1.67 2xm

0.167 0.333 o*soo 0.667 0.833 1.00

0.167 0.333 0.500 0.667 0.833 1.00

o-0 o-0 0.0 0.0 0.0 0.0

D

0.333 0.667 1.00 1.33 1.67 2.00

0.145 0.293 0.459 0.655 0.828 1.00

0.167 0.333 0.500 0.667 0.833 1.00

-13.2 -12.0 -8.2 -1.8 -0.6 0.0

E

0,333 0.667 1.00

0.333 0.667 1.00

0.333 0.667 1.00

0.0 0.0 0.0

A

Radius

1.00

B

1.00 1.33

% Error

436

R. W. STEPHENSON

Figure

5. Finite element

ET AL.

model of Figure 2 configuration.

TABLE 3 Axial displacements of beam in Figure 5 Axial location

Non-dimensional displacement

Radius

Beam theory displacement

% Error

A

0.333 0.677 1.00

0.333 0.667 1.00

0.333 0.667 1.00

0.0 o-0 0.0

B

0.333 0.667 1.00

0.333 0.677 1.00

0.333 0.667 1.00

0.0 0.0 0.0

C

0.333 0.667 1.00 1.33 1.67 2.00

0.479 0.908 1.08 0.826 0.904 1.00

0.167 0.333 0.500 0.667 0.833 1.00

187.0 173.0 116.0 23.8 8.5 0.0

D

0.333 0.667 1.00 1.33 1.67 2.00

0.175 0.343 0.507 0.668 0,832 1.00

0.167 0.333 0.500 0.667 0.833 1.00

4.8 3.0 1.4 0.1 -1.1 0.0

E

0.333 0.667 1.00 1.33 1.67 1.00

0.166 0.333 0.499 0.666 0.833 1.00

0.167 0.333 0.500 0.667 0.833 1.00

0.6 0.0 0.2 0.1 0.0 0.0

at the small diameter (free) end. The large diameter end is fixed. For this geometry, the deviation from beam theory is quite significant at the diameter change. The distortion of the cross-section is greatest at the location of the diameter change, decreases along the larger diameter, and becomes essentially zero after about one diameter of length. 3.2.

DYNAMIC

ANALYSIS

OF

A

ROTOR

[18] used an experimental three-disk rotor to compare predicted and measured frequencies for a free-free (non-rotating) situation. Table 4 gives a summary of the first five free-free natural frequencies obtained by using (a) Vance’s experimental data, (b) Vance’s transfer matrix calculations, along with calculations of this paper obtained by using (c) beam finite elements, (d) axisymmetric solid elements from ANSYS, and (e) the cubic axisymmetric elements developed in reference [14]. Vance,

Murphy

and

Tripp

ROTOR MODELLINGWITH SOLID

437

ELEMENTS

TABLE 4 Summary of free-free natural frequencies of three-disk laboratory rotor [18], elastic modulus = 2.068 E 11 N/m*, mass density = 20.27 kg/m’, Poisson ratio = O-25, percent error in parentheses Mode , Measured (Hz) Transfer matrix, lumped mass: (1) with shear, without correction factors (2) with shear, with correction factors

Finite element, beams with consistent mass: (1) with shear, without correction factors (2) without shear, without correction factors

1

2

3

4

\ 5

94

201

356

463

832

95 (1.1)

207 (0.0)

353 (-0.8)

408 (-12)

850 (2.2)

200 (-3.5)

345 (-3.1)

403 (-13)

786 (-5.5)

101 (8.0)

223 (7.7)

383 (7.6)

491 (6.0)

973 (17)

102 (8.5)

225 (8.7)

389 (9.3)

494 (6.7)

1117 (34)

(3q77)

214 (3.5)

361 (3.1)

488 (5.4)

890 (7.0)

206 (-0.5)

355 (-0.3)

460 (-0.6)

839 (0.8)

208 (0.5)

355 (-0.3)

471 (1.7)

931 (12)

91 (-2.9)

(3) with shear, with correction factors

Finite element, axisymmetric from ANSYS

solid elements

Finite element, cubic axisymmetric

elements

93 (-1.1) 95 (1.1)

3.2.1. Beam results The transfer matrix results of reference [ 181, obtained by using a lumped mass formulation, appear to give reasonably good agreement, at least for the first three modes. However, the model of reference [ 181 does not include any correction factors for the section change. If the correction factors are included (as they should be), the frequencies are significantly lower than measured. The authors believe that the good agreement without the correction factors is due to two compensating effects: the lumped mass formulation in this case gives frequencies that are too low, while neglecting the correction factors increases the stiffness and the calculated frequencies. The finite element results with consistent mass beam elements show that the error from ignoring the approximate correction factor is greater than the error due to ignoring shear. 3.2.2. Solid model results The first model discussed is based on the ANSYS STIF42 element. This is a four-node isoparametric element, with extra shape functions added for improved bending behavior. This element has a consistent mass matrix, but is non-conforming because of the extra shape functions used. (Later versions of the ANSYS program have a modified extra shape function, which passes the patch test.) The element provides a displacement function and performance approaching that of a conforming quadratic element, because of the extra shape functions employed. The model of Figure 6 is fairly detailed, and should provide a model of bending behavior of the rotor which is very close to the physical structure.

438

R. W. STEPHENSON

Figure 6. Finite element model of three-disk

ET AL.

laboratory

rotor with ANSYS

elements.

(Note that even this model will ignore the very local strain distribution in the sharp corners at the diameter steps; if these were of interest, the level of modelling detail could be further extended to include these effects. For our purposes, these local concentrations are of much less importance than the distortions adjacent to the corners.) The natural frequencies of the ANSYS model were evaluated for free boundary conditions, which were approximated in the test by Vance et al. [18]. This gives two zero-frequency rigid body modes. The beam-like form of the rotor leads to bending modes which approximate the typical first, second and third modes of beam theory, with variations due to the large disks in the center of the rotor. The matrix reduction method was used in the ANSYS solution, with 23 degrees of freedom being retained. These master degrees of freedom were distributed along the length of the rotor at the rotor centerline. Vance reported his results to the nearest hertz, so that round-off could lead to some variation in the error comparison of Table 4. However, the first three modes of the ANSYS model are within 1 Hz of the experimental results, and the percent error is within about 1% for the first five modes. The first four modes are slightly lower than experimental, but it is difficult to judge if this is due to some physical restraint in the experimental set-up, or to convergence characteristics in the finite element model. The agreement of the ANSYS model in Table 4 is considered to be an excellent corroboration of both the finite element method and the experimental techniques of Vance, Murphy and Tripp. The model was also prepared with a cubic element formulation in a code developed by the authors. The element has the same basic Fourier series approach as the ANSYS element, but is conforming and has a subparametric cubic formulation. The four corner nodes define the element geometry, allowing only straight edges for the model discretization. This should not be a limitation for the modelling of most rotors. To define the element deflection, the element requires 12 nodes, with two mid-side nodes on each edge so that each edge can then displace in a cubic shape. The higher order should provide greater accuracy for a given number of elements. Figure 7 shows the finite element model used for the cubic elements. In preparing this model, the absolute minimum number of elements was used which would still represent the geometry of the rotor. Since only four-sided elements were developed, the model consists of a single layer of elements for the basic shaft portion, with a second layer for radial projections on the shaft such as the three disks. Master degrees of freedom were selected in a manner similar to those in the ANSYS model, with a total of 25 retained. The fairly coarse mode1 with cubic elements performs nearly as well as the detailed ANSYS model for the first four modes. The fifth mode is in greater error (12%), apparently because the cubic displacement in the axial direction is not adequate for the bending mode shape in the higher modes. It is also possible that the cubic shape in the radial direction provides a lower degree of accuracy than the larger number of ANSYS elements used in the first model.

Figure 7. Finite element model of three-disk

laboratory

rotor with cubic axisymmetric

elements.

ROTOR

MODELLING

WITH

SOLID

ELEMENTS

439

4. CONCLUSIONS

1. The distortion of the cross-section of a shaft that occurs at a stepped diameter change tends to decrease the shaft bending stiffness and natural frequencies calculated from classical beam theory. 2. The results illustrate the need for using some type of a correction factor when beam elements are used to model shafts with significant diameter changes. 3. The error that occurs from ignoring the section change will depend on the severity of the diameter change. For the cases studied here, this error is greater than the error that occurs from neglecting shear deflection. 4. A review of the results in reference [18], which reported relatively good agreement between lumped mass and transfer matrix results, indicates that the lumped mass approximation is actually lowering the frequencies, somewhat compensating for the increase due to the “plane section” assumption. 5. The excellent agreement between the measured frequencies and the frequencies calculated by using an axisymmetric solid finite element model shows the even greater accuracy possible with these elements. The dynamic reduction procedure gives excellent results when only a small fraction of the total degrees of freedom are retained. 6. A relatively coarse model with the axisymmetric cubic elemen: (described in the Appendix) can be used to give acceptable accuracy as well. Extension of the formulation to include gyroscopic effects is necessary for the rotating case-this work is in progress. ACKNOWLEDGMENT

This research was partially supported with funds from the National Science Foundation (Grant No. RII-8610671) and the Commonwealth of Kentucky through the Kentucky EPSCoR Program. REFERENCES

and F. K. ORCUI-~ 1967 American Society of Mechanical Engineers Transactions, Journal of Engineering for Industry 89, 785-796. Calculation and experiments on the unbalance

1. J. W. LUND

response of a flexible rotor. 2. R. UHRIG 1966 Journal of Sound and Vibration 4(2), 136-148. The transfer matrix method seen as one method of structural analysis among others. 3. R. RUHL 1970 Ph.D. TIesis, Cornell University. Dynamics of distributed parameter rotor systems: transfer matrix and finite element techniques. 4. K. E. ROUCH and J. S. KAO 1979 Journal of Sound and Vibration 66(l), 119-140. A.tapered beam finite element for rotor dynamics analysis. 5. L.M. GREENHILL,W.B.BICKFORD~~~ H.D.NELSON 1985 AmericanSocietyofMechanical Engineers Transactions, Journal of Vibration, Acoustics, Stress, and Reliability in Design 107, 421-427. A conical beam finite element for rotor dynamics analysis. 6. K. E. ROUCH and J. S. KAO 1980 American Society of Mechanical Engineers Transactions, Journal of Mechanical Design 102,360-368. Dynamic reduction of rotor dynamics by the finite

element method. 7. R. T. SEVERN 1970 Journal of Strain Analysis 5, 239-241. Inclusion of shear deformation in the stiffness matrix for a beam element. 8. H. D. NELSON and J. M. MCVAUGH 1976 American Society of Mechanical Engineers Transactions, Journal of Engineering for Industry 98,593-600. Dynamics of rotor-bearing systems using finite elements. 9. R. L. ESHELMAN and R. A. EUBANKS 1969 American Society of Mechanical Engineers Transactions, Journal of Engineering for Industry 91, 1180- 1188. On the critical speeds of a continuous rotor. 10. D. L. THOMAS, J. M. WILSON and R. R. WILSON 1973 Journal of Sound and Vibration 31, 315-330. Timoshenko beam finite elements.

440

R. W. STEPHENSON

ET

AL.

11. M. J. BROWN 1952 Ph.D. Thesis, Purdue University. Lateral and torsional stiffness of shafts having circular cross section and varying diameters. 12. R. D. COOK 1974 Concepts and Applications of Finite Element Analysis, New York: John Wiley. 13. M. GERADIN and N. KILL 1984 Engineering Computations 1, 52-64. A new approach to finite 14.

element modelling of flexible rotors. R. ARORA 1987 M.S. Thesis, University of Kentucky. Axisymmetric finite elements for rotor dynamics.

15. G. J. DESALVO and J. A. SWANSON 1987 ANSYS

Users Manual. Swanson Analysis Systems, Inc. 16. R. J. GUYAN 1965 American Institute of Aeronautics and Astronautics Journal 3(2), 380. Reduction of stiffness and mass matrices. 17. R. W. STEPHENSON and K. E. ROUCH 1988 Second International Symposium on Transport Phenomena, Dynamics, and Design of Rotating Machinery, 3-6 April, 1988, Honolulu, Hawaii.

Rotor modeling considerations with shaft diameter changes. 18. J. M. VANCE, B. T. MURPHY and H. A. TRIPP 1985 ASME Paper VIB85RD-142, presented at the Design Engineering Division Conference, lo-13 September, 1985, Cincinnati, Ohio. Critical

speeds of rotating machinery: computer predictions The rotor mass-elastic model.

vs. experimental

measurements,

Part I.

APPENDIX: MASS AND STIFFNESS MATRICES OF AN AXISYMMETRIC ELEMENT The development is briefly given here for a subparametric cubic quadrilateral axisymmetric element, shown in Figure 8. Four nodes are used to define the geometry shape function, and 12 nodes are used to define the displacements.

Figure 8. Configuration of subparametric cubic axisymmetric element.

Lagrangian interpolation functions are used along the 5 and n directions geometry. The global radial and axial co-ordinates, r and z, are expressed as i=4

i=4 r =

for the

C N,,ri, i=l

z=C i=l

Nz. & ‘3

(AlI

where Ns, are the geometry shape functions Ns, = I/4(1 - 5)(I - n),

Ngz= l/4(1 + &)(I - n),

N,,=I/4(1+5)(l+n),

N&=1/4(1--&)(l+n).

(A2)

At each node there is a radial, axial, and circumferential degree of freedom. These displacements are approximated from all 12 nodes and corresponding displacement shape functions, Ni, as i=12

24;= C NiUz COS(no), i=l

i=l2

u; =

1 i=l

i=I2

Niti:, cos (no),

u:=

1 - Niuz;sin (tie). i=l

(A3)

ROTOR

MODELLING

WITH

SOLID

441

ELEMENTS

For lateral bending of a beam, the only mode of interest is n = 1. Equation (A3) reduces to isI2

i=I2

i=12

U, =

1 NiU,, COS0,

U, =

i=l

I.&= 1 - NiU,, sin 8.

1 NiU, COS0, i=l

(A4)

i=l

Equation (A4) can be written in matrix notation as N, cos &J

0 0

0

0

N, cos 8 0

-N, sin 8

0

(A5) a.*

or {u} = [N](A). The displacement shape functions are derived in a straightforward corner nodes 1,2,3,4, where .$=*l and n=fl,

manner. For the

Ni=(1/32)(1+~~i)(l+nni)[-10+9(~2+~2)]. For the vertical midside nodes 7,8,11,12,

(A6a)

where 5 = *l and n = *l/3, (A6b)

Ni = (9/32)(1+ 551)(1- n2)(l +9TTi)* For the horizontal midside nodes 5,6,9,10,

where 5 = *l/3

and 7 = *l, (A6c)

Ni=(9/32)(l+nTi)(l-q2)(1+955i)*

Since non-dimensionalized co-ordinates are used to define shape functions, it is necessary to express the global derivatives in terms of the local derivatives as (A7) where the Jacobian [J] is given by

[Jl= From equation (Al) and differentiating -(l-n) -(l-t)

az/@

ar/&$

az/a7)

ar/ag

(A81

I ’

(A2) the Jacobian can be shown to be

(1-n) -(l+[)

(l+~) (l+t)

-(l+n) (1-t)

(A9)

The local derivatives are defined as follows. For nodes 1,2,3 and 4, alv,la5=(1l32)(1+~~i)[185(1+~~i)+~i(-10+9~2+P-qf~_ (AlOa)

aNi/a~=(1/32)(1+&i)[18n(l+~ni)+~i(-10+9~2+n2)]. For nodes 5,6,9 and 10, aNilX=

(g/32)(1

i3Nilaq

+ TTi)(95i-25-27t25i),

=(9/32)(1

-t2)(1+9&)+ (AlOb)

For nodes 7,8,11 and 12, aNilay=(9/32)(1-712)(1+97717i),

6’Nilt3q

= (9/32)(1+ Tni)(9ni -27) -27~7~7~). (AlOc)

442

R.

W. STEPHENSON

ET AL.

To obtain the element stiffness matrix, consider the strain vector {E}, alFIr 0

0

0

d/dZ

0

l/r a/az

alar

0

l/r a/at9

0

l/t-a/at3

[Ll{uI,

=

0

(All)

alar-l/r

and the stress vector, {a} = [C](E), where [C] is the constitutive matrix which is a function of the elastic modulus E and the Poisson ratio V. The stiffness matrix for an element can then be shown to be [K]=j

[BITIC][B]dV=

de dr dz,

[BITIC][B]r

(Al21

V

where [B]=[L][N] (aN,/ar)

cos

0

8

0 (N,/r)

=

(aN,/az) cos

(aN,/az)

8

0

e

f.. similar ... for . . . nodes *** 2to1;

0

sin 0

...

cos

8

0

sin 8

-( N,/r)

sin 8

N,/r)

cos

...

-( N,/r)

(aNJar)

e

0 -(

e

0

8

cos

cos

0

-(aN,/az) -(aN,/ar

cos

- N,/r)

All entries of the matrix product

[BIT[ C][ B] contain either co? 8 or sin2 8 terms. Integration can be carried out with respect to 8 independently of the other co-ordinates and upon noting that Ii” co? 0 de =,irr sin* 8 de = V, the resulting stiffness matrix is of the form

EkJMr dr dz

[K]=[&r/(l+v)(l-2v)]

where aN,/ar

0

0

[k&x6 =

aN,/az

0

0

NJ r

aN,/az

0

aNJar

-NJr

0

similar for nodes 2 to 12 and I- ,

(1 -v)

aN,/ar+

v(aN,/ar+

vN,/r N,/r)

~dN,/ar+(l-v)N,/r (i-2~)/2 aN,/az 0

(1-2v)/2

N,/r

Y aN,/az (I-

V) aN,/az

Y aN,/az (1-2v)/2 aNJar (1-2v)/2 N,/r 0

0 N,/ r aN,/az

(A13)

ROTOR

MODELLING

WITH

SOLID

443

ELEMENTS

-vN,/r

I*. Similar-

-vN,/r

...

-(1-v)N,/r

..e *. .

0 (1-2v)/2

fW,/az

(l-2v)/2(aN,/ar-

N,/r)

for nodes

**.

2 to

.a.

12

.

_

The integration over the element surface is performed by using the mapping dr dz = det J dt dv with integration limits -1 to +l. The element mass matrix is obtained in the normal manner from [M]=[

p[NITIN]dV= V

111

p[ NIT[ N]r de dr dz,

(A14)

where p is the mass density and [N] is given in equation (A5). Integrating with respect to 6 as beforegives [M]=p~jj[m]rdrdz, where N:

0

0

N:

0

0

N,;N,

0

[ml = 0 0

N,N,2

0

0

,..

0

N,N,z

.. . . ..

0

0

0 N,N,,

iI

0

a..

q-. N,;N,*

NzN

...

0

0

. ..

0

N,JG 0

’

0 N,&,J