Journal of Sound and Vibration (1995) 187(5), 841–849
VIBRATIONS OF TIMOSHENKO BEAMS BY VARIABLE ORDER FINITE ELEMENTS A. H Institute of Mechanical Engineering, University of Tlemcen, Tlemcen 13000, Algeria (Received 28 October, and in final form 17 January 1995) This paper presents a four-node Timoshenko beam finite element with variable degrees of freedom. Both the element transverse displacement and the rotation of the beam cross-section are described by a cubic polynomial plus a variable number of trigonometric sine terms. The polynomial terms are used to describe the transverse displacements and the rotations of the beam cross-section at the element’s four nodes and the sine terms are used to provide additional freedom to the interior of the element. The four nodal transverse displacements and rotations of the beam cross-section and the amplitudes of the trigonometric sine terms are used as generalized co-ordinates. Inter-element compatibility is achieved by matching the generalized co-ordinates at the element end nodes. Numerical results of frequency calculations are given for simply supported and cantilever beams with two different slenderness ratios. Comparisons are made with exact Timoshenko beam solutions and with finite element solutions for the degenerate case with no trigonometric terms to represent a polynomial finite element. Comparisons show that using one or two variable order Timoshenko beam elements with a few trigonometric terms yields a better accuracy with fewer system degrees of freedom than using many polynomial Timonshenko beam finite elements. 71995 Academic Press Limited
1. INTRODUCTION
This paper deals with a finite element with variable degrees of freedom for the vibration of Timoshenko beams. Milsted and Hutchinson [1] and Hutchinson and Benitou [2] have shown that elements based on similar formulations for membrane and plate vibrations were very successful, and that provided the motivation for this work. Timoshenko theory is well-known for flexural vibrations of beams. This theory takes into account the effects of two major factors, namely rotary inertia and shear deformation. The theory which takes into account the effect of rotary inertia alone was first developed by Rayleigh [3]. Timoshenko [4, 5] extended this theory to include shear deformation. Exact solutions based on the Timoshenko theory are available for beams with simply supported ends [6] and with other boundary conditions [7]. Comparisons of the Timoshenko theory with the elasticity theory, as well as useful suggestions of shear coefficients that give good matching with the elasticity solutions, are available for beams of circular cross-section [8] and for beams of rectangular cross-section [9]. A number of Timoshenko beam finite elements in which only polynomial terms are used to describe the element degrees of freedom are available in the literature. Dawe [10] presented a new Timoshenko beam finite element and reviewed the other existing elements. The idea of using trigonometric terms in the finite element method is not new. Pian [11] described the concept of using more co-ordinates than the element nodal displacements in 841 0022–460X/95/450841+09 $12.00/0
7 1995 Academic Press Limited
842
.
deriving element stiffness matrices. Krahula and Polhemus [12] used the Fourier series for a rectangular plane stress element. The desirability of using high order finite elements for vibration problems has been shown by many authors. Thomas and Dokumaci [13] have shown that high order finite elements yield improved results for the vibration of tapered beams. On the basis of a work on the vibration of plates using Mindlin theory which is equivalent in plate analysis to that of Timoshenko theory in beam analysis, Dawe [14] suggested that an increase in efficiency would result if finite elements of order higher than those used in previous Timoshenko beam models were utilized. In this paper a four-node finite element with variable degrees of freedom is developed for obtaining natural frequencies and mode shapes of Timoshenko beams. Both the element transverse displacement w and the rotation of the beam cross-section u are described by a cubic polynomial plus a variable number of trigonometric sine terms. The polynomial terms are used to describe the degrees of freedom at the element’s four nodes located at each end and at the one-third points and the sine terms are used to provide additional freedom to the interior of the element. If the polynomial terms are used by themselves they will be sufficient to describe the displacement and the rotation of the beam cross-section within the element and therefore will be a good choice for a standard type of finite element. The polynomial terms used in this paper lead to a Timoshenko beam element similar to the one developed by Carnegie et al. [15]. The additional trigonometric sine terms allow for a better description of the transverse displacement and the rotation of beam cross-section within the element. The transverse displacements and rotations of the beam cross-section at the four nodes and the amplitudes of the trigonometric terms are used as generalized co-ordinates. Inter-element compatibility is achieved by matching the generalized co-ordinates at the element end nodes. The frequencies for the four lowest modes for two simply supported prismatic beams and the frequencies for the three lowest modes for two cantilever prismatic beams are found. These particular examples are chosen because exact solutions based on Timoshenko theory are available for comparisons. Comparisons are also made with finite element solutions for the degenerate case with no trigonometric terms to represent the element, which will be referred to as the polynomial element to distinguish it from the trigonometric element. The formulation in this paper is not limited to isotropic linearly elastic prismatic beams, but could be easily modified to account for the variation of the depth of the beam cross-section or the variation of the material properties. 2. FORMULATION
The co-ordinate system used to define the geometry of a four-node Timoshenko beam element is shown in Figure 1. The local and global co-ordinates are related by the equation (a list of notation is given in the Appendix) j=x/L.
Figure 1. Element co-ordinates.
(1)
843
Figure 2. Deformations of the beam differential element.
The potential energy PE and the kinetic energy KE of the prismatic Timoshenko beam element are PE=
$ g0 1
1 EI 2
1
0
KE=
2
1 du L dj+kGA L dj
$g
rv 2 A 2
g0 1
0
g
1
0
1 dw L dj , L dj
%
1
w 2L dj+I
1 % 2
u−
u 2 L dj .
0
(2)
(3)
Generic displacements for this element are the transverse displacement w and the rotation of beam cross-section u as shown in Figure 2. The choice of u as a second degree of freedom instead of c or dw/dx allows a better representation of all possible beam boundary conditions as well as geometric discontinuities. The displacement functions assumed for this element are w=q1+q3 j+q5 j 2+q7 j 3+qm1 sin am j,
(4)
u=q2+q4 j+q6 j 2+q8 j 3+qm2 sin am j,
(5)
where am=mp,
m=1, 2, . . . , M,
(6)
m2=m1+1.
(7, 8)
and the indices are defined as m1=2m+7,
The chosen polynomial displacement functions are complete and balanced and represent an obvious choice for a beam element. The polynomial terms in the assumed displacement field are used to define the transverse displacements and the rotations of beam cross-section at the element’s four nodes and the trigonometric sine terms are used to give additional freedom to the interior of the element. The quantities needed to form the element stiffness and mass matrices are obtained in matrix form and are w=Cq,
u=Dq,
(1/L) du/dj=Fq,
u−(1/L) dw/dj=Hq,
(9–12)
where C=[1, 0, j, 0, j2, 0, j3, 0, sin am j, 0],
(13)
D=[0, 1, 0, j, 0, j 2, 0, j 3, 0, sin am j],
(14)
.
844
F=[0, 0, 0, 1/L, 0, 2j/L, 0, 3j 2/L, 0, (am /L) cos am j],
(15)
H=[0, 1, −1/L, j, −2j/L, j 2, −3j 2/L, j 3, −(am /L) cos am j, sin am j],
(16)
qT=[q1 , q2 , q3 , q4 , q5 , q6 , q7 , q8 , qm1 , qm2 ].
(17)
and
Substituting equations (9)–(12) into equations (2) and (3) gives KE=12 v 2qTMq q,
PE=12 qTKq q,
(18, 19)
where Kq and Mq are, respectively, the element stiffness and mass matrices, expressed in the q co-ordinate system as Kq=EI
g
1
0
FTFL dj+kGA
g
1
HTHL dj,
Mq=rA
0
g
1
CTCL dj+rI
0
g
1
DTDL dj.
0
(20, 21) The element stiffness matrix Kq and mass matrix Mq expressed in the q co-ordinate system will be obtained by integrating explicitly the expressions given by equations (20) and (21). A new set of generalized co-ordinates p is chosen in order to satisfy inter-element compatibility. The relation between the p co-ordinates and the q co-ordinates is found by applying the element ‘‘boundary conditions’’ and is q=Tp,
(22)
pT=[w1 , u1 , w2 , u2 , w3 , u3 , w4 , u4 , wm , um ].
(23)
where
The generalized co-ordinates w1 , u1 , w2 , u2 , w3 , u3 , w4 , and u4 are the transverse displacements w and the rotations of the beam cross-section u at the element’s four nodes and the generalized co-ordinates wm and um are the amplitudes of the trigonometric functions for the transverse displacement w and the rotation of the beam cross-section u in the interior of the element. These generalized co-ordinates are shown in Figure 3. Because the polynomial terms are used to describe the transverse displacements and the rotations of the beam cross-section at the element’s four nodes and the amplitudes of the trigonometric sine terms are used to provide additional freedom for the transverse displacement and the rotation of the beam cross-section in the interior of the element the transformation from the q co-ordinate system to the p co-ordinate system remains uncoupled regardless of the number of trigonometric terms used. The transformation matrix T is
Figure 3. Element generalized co-ordinates.
845
constituted of a fixed diagonal 8×8 block and a variable number of diagonal 2×2 blocks and zero coefficients outside these blocks. The matrix T takes the following form for M=1:
K 1 0 0 0 0 0 0 0 L G 0 G 1 0 0 0 0 0 0 G −11/2 G 0 1 0 9 0 −9/2 0 G G −11/2 0 1 0 9 0 −9/2 G 0 G G 9 0 −9/2 0 −45/2 0 18 0 G . T= G 0 9 0 −9/2 0 −45/2 0 18 G G −9/2 G 0 9/2 0 27/2 0 −27/2 0 G 0 G −9/2 0 9/2 0 27/2 0 −27/2 G G 1 0G G k 0 1l (24) The order N of the transformation matrix T depends on the number of trigonometric terms M chosen. The matrix T is found by simply expanding it along the diagonal with as many 2×2 blocks as the number of terms used. The order N of the element stiffness, mass, and transformation matrices is (25)
N=2M+8.
The element stiffness and mass matrices are transformed to the p co-ordinate system by using the relations Kp=TTKq T,
Mp=TTMq T.
(26, 27)
The element stiffness matrix Kp and mass matrix Mp are then assembled into the system stiffness and mass matrices by direct summation. Any of the known techniques that solves a generalized eigenvalue problem can then be used to find the frequencies (eigenvalues) and mode shapes (eigenvectors). 3. NUMERICAL RESULTS
Results of the application of the variable order Timoshenko beam element to the calculation of the frequency parameters V for the four lowest modes of two simply supported beams and the three lowest modes of two cantilever beams are given respectively in Tables 1
T 1. Frequency parameters V for the four lowest modes of two simply supported beams Exact L/r Mode solution
1T2 (10)
1T3 (12)
1T4 (14)
1T5 (16)
2T2 (20)
2T3 (24)
5P (30)
6P (36)
7P (42)
12·5
1 2 3 4
0·09122 0·29371 0·53144 0·77774
0·09122 0·29424 0·58058 0·93542
0·09122 0·29372 0·53328 0·86756
0·09122 0·29371 0·53147 0·78125
0·09122 0·29371 0·53146 0·77781
0·09122 0·29371 0·53151 0·77846
0·09122 0·29371 0·53145 0·77775
0·09122 0·29373 0·53165 0·77902
0·09122 0·29372 0·53155 0·77819
0·09122 0·29371 0·53145 0·77793
25
1 2 3 4
0·02469 0·09122 0·18487 0·29371
0·02469 0·09149 0·21409 0·39596
0·02469 0·09122 0·18606 0·34575
0·02469 0·09122 0·18489 0·29644
0·02469 0·09122 0·18488 0·29377
0·02469 0·09122 0·18491 0·29423
0·02469 0·09122 0·18487 0·29372
0·02469 0·09123 0·18500 0·29454
0·02469 0·09122 0·18492 0·29401
0·02469 0·09122 0·18488 0·29384
.
846
T 2. Frequency parameters V for the three lowest modes of two cantilever beams Exact L/r Mode solution
1T2 (10)
1T3 (12)
1T4 (14)
1T5 (16)
2T2 (20)
2T3 (24)
5P (30)
6P (36)
7P (42)
20
1 2 3
0·01378 0·01378 0·01378 0·01378 0·01378 0·01378 0·01378 0·01378 0·01378 0·01378 0·07503 0·07506 0·07504 0·07503 0·07503 0·07503 0·07503 0·07503 0·07503 0·07503 0·17988 0·18175 0·18002 0·17989 0·17988 0·17988 0·17988 0·17991 0·17990 0·17989
50
1 2 3
0·00226 0·00226 0·00226 0·00226 0·00226 0·00226 0·00226 0·00226 0·00226 0·00226 0·01378 0·01380 0·01378 0·01378 0·01378 0·01378 0·01378 0·01378 0·01378 0·01378 0·03712 0·03794 0·03723 0·03712 0·03712 0·03712 0·03712 0·03714 0·03713 0·03712
and 2. The frequency parameter V used herein is dimensionless and is obtained by multiplying the natural frequency v by r and dividing it by the shear wave velocity. In Tables 1 and 2 the columns labelled nTm refer to solutions with n variable order Timoshenko beam elements in the complete beam with m trigonometric terms in each element. The columns labelled nP refer to solutions with n Timoshenko beam elements with no trigonometric terms in the complete beam. The numbers in parenthesis refer to the numbers of system degrees of freedom excluding the restrained ones. Because the beam has no geometric or material discontinuities it is first discretized into one element and the number of trigonometric terms M in each element is varied. The beam is also discretized into two elements to show that more than one single element can be used in the complete beam. The results obtained are for n=0·3 for all beams, k=0·85 for the simply supported beams, and k=0·65 and 2/3 for L/r=20 and 50 respectively for the cantilever beams. The results of Tables 1 and 2 show that the variable order finite element solutions always converge to the exact solutions. In order to have a good idea of the manner of convergence of the one-element solution the frequency parameter is plotted versus the number of trigonometric terms in Figures 4 and 5 for the simply supported beams and in Figures 6 and 7 for the cantilever ones. These figures confirm that the variable order finite element solutions converge from above to the exact solutions represented by dashed lines. As can be seen, convergence to the exact values is in fact very rapid for all beams even for the very deep ones in which shear deformation is more significant. As can be seen convergence to the exact solutions in the lowest two modes is much faster than in the third or fourth modes. This is due to the fact that the higher is the mode the more complex is the mode shape. Therefore the higher the mode is, the more trigonometric terms are needed to describe accurately this mode. It can also be seen that for the cantilever beams only two trigonometric terms were
Figure 4. Convergence to V1, V2 , V3 , and V4 , of the simply supported beam with L/r=12·5.
847
Figure 5. Convergence to V1, V2 , V3 , and V4 , of the simply supported beam with L/r=25.
sufficient to obtain adequate accuracy for the third mode. However, for the simply supported beams three trigonometric terms were needed to obtain adequate accuracy for the third mode. The fact that the cantilever beam has only one restrained end makes its mode shapes simpler than those of the simply supported one. Therefore for the cantilever beams fewer trigonometric terms are needed to describe accurately a particular mode. Comparisons of the results of Tables 1 and 2 show that by using one or two variable order Timoshenko beam elements in the complete beam and varying the number of trigonometric terms highly accurate answers are obtained regardless of the beam’s slenderness ratio. It is also shown in Tables 1 and 2 that the two-element solution requires more system degrees of freedom than the one-element solution to obtain equivalent accuracy. The two-element solution is, however, computationaly more efficient than the one-element solution. In fact, in the two-element solution the global stiffness and mass matrices have smaller bandwidths than those for the one-element solution and therefore require less storage space. Furthermore, it is also shown in Tables 1 and 2 that the variable order Timoshenko beam element produces a better accuracy than the polynomial element, with fewer system degrees of freedom. The computer time required for the execution of the transformations in equations (26) and (27) can be greatly reduced if the zeros outside the diagonal blocks of the matrix T are disregarded during multiplication.
Figure 6. Convergence to V1, V2 , and V3 , of the cantilever beam with L/r=20.
848
.
Figure 7. Convergence to V1 , V2 , and V3 of the cantilever beam with L/r=50.
4. DISCUSSION AND CONCLUSIONS
The variable order Timoshenko beam finite element presented in this paper is formulated in terms of a cubic polynomial and a variable number of trigonometric sine terms for both the transverse displacement and the rotation of the beam cross-section. One advantage of this element is that highly accurate frequencies for Timoshenko beams can be obtained with a small number of system degrees of freedom. In comparison to a Timoshenko beam finite element with no trigonometric terms, it showed a better accuracy with fewer system degrees of freedom. The element formulation can be slightly modified to account for variations of the depth of the beam cross-section or the material properties. One possible drawback of this element is that it requires extra matrix multiplications related to the transformation of system coordinates. However, the extra computer time associated with the transformation procedure can be greatly reduced if zeros outside the diagonal blocks of the transformation matrices are disregarded during multiplication. A major advantage of this element is that one can vary the number of system degrees of freedom in a structure without changing the gridwork of elements. Results can then be obtained to any desired degree of accuracy by simply increasing the numbers of trigonometric terms in each element. The polynomial element does have the advantage that the transformation of co-ordinates is not required. While the variable order Timoshenko beam element does not have this advantage, the higher accuracy obtained by using a smaller number of system degrees of freedom more than makes up for this deficiency. REFERENCES 1. M. G. M and J. R. H 1974 Journal of Sound and Vibration 32, 327–346. Use of trigonometric terms in the finite element method with application to vibrating membranes. 2. J. R. H and J. J. B 1977 Journal of the Engineering Mechanics Division, American Society of Civil Engineers 103, 779–787. Variable order finite elements for plate vibration. 3. L R 1877 Theory of Sound (two volumes). New York: Dover; second edition, 1945 re-issue. 4. S. P. T 1921 Philosophical Magazine 41, 744–746. On the correction for shear of the differential equation for transverse vibrations of prismatic beams. 5. S. P. T 1922 Philosophical Magazine 43, 125–131. On the transverse vibration of bars of uniform cross-section. 6. C. D 1954 Quarterly of Applied Mathematics 12, 175–187. On the Timoshenko theory of transverse beam vibrations.
849
7. T. C H and C. S. K 1962 Developments in Theoretical and Applied Mechanics 1, 59–71. New tables of eigenfunctions representing normal modes of vibration of Timoshenko beams. 8. J. R. H 1981 Journal of Applied Mechanics 48, 923–928. Transverse vibrations of beams, exact versus approximate solutions. 9. J. R. H and S. D. Z 1986 Journal of Applied Mechanics 53, 39–44. On the transverse vibration of beams of rectangular cross-section. 10. D. J. D 1978 Journal of Sound and Vibration 60, 11–20. A finite element for the vibration analysis of Timoshenko beams. 11. T. H. H. P 1964 American Institute of Aeronautics and Astronautics Journal 2, 576–577. Derivation of element stiffness matrices. 12. J. L. K and J. F. P 1968 American Institute of Aeronautics and Astronautics Journal 6, 726–728. Use of Fourier series in the finite element method. 13. J. T and E. D 1973 Aeronautical Quarterly 24, 39–46. Improved finite elements for vibration analysis of tapered beams. 14. D. J. D 1978 Journal of Sound and Vibration 59, 441–452. Finite strip models for vibration of Mindlin plates. 15. W. C, J. T and E. D 1969 Aeronautical Quarterly 20, 331–332. An improved method of matrix displacement analysis in vibration problems. APPENDIX: NOTATION E G r n k A I L r x j w c u
modulus of elasticity shear modulus mass density Poisson ratio shear coefficient beam cross-sectional area beam moment of inertia beam length =zI/A , radius of gyration of beam cross-section element global co-ordinate element dimensionless co-ordinate, origin at node 1 transverse displacement shear distorsion =c+dw/dx, rotation of beam crosssection
element stiffness matrix in the q coordinate system Mq element mass matrix in the q co-ordinate system Kp element stiffness matrix in the p coordinate system Mp element mass matrix in the p co-ordinate system T transformation matrix M number of trigonometric terms PE potential energy KE kinetic energy v natural frequency V =vr/zG/r , non-dimensional frequency
Kq