- Email: [email protected]
A new method to determine the shear coefficient of Timoshenko beam theory
Journal of Sound and Vibration 330 (2011) 3488–3497
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
A new method to determine the shear coefficient of Timoshenko beam theory K.T. Chan a, K.F. Lai b, N.G. Stephen c, K. Young b, a b c
Department of Mechanical Engineering, The Hong Kong Polytechnic University, Hong Kong, PR China Department of Physics, The Chinese University of Hong Kong, Hong Kong, PR China School of Engineering Sciences, Mechanical Engineering, University of Southampton, Highfield, Southampton SO17 1BJ, UK
a r t i c l e in f o
abstract
Article history: Received 17 November 2010 Received in revised form 8 February 2011 Accepted 8 February 2011 Handling Editor: S. Ilanko Available online 5 March 2011
The frequencies o of flexural vibrations in a uniform beam of arbitrary cross-section and length L are analysed by expanding the exact elastodynamics equations in powers of the wavenumber q ¼ mp=L, where m is the mode number: o2 ¼ A4 q4 þ A6 q6 þ . The coefficients A4 and A6 are obtained without further assumptions; the former captures Euler–Bernoulli theory while the latter, when compared with Timoshenko beam theory rendered into the same form, unambiguously yields the shear coefficient k for any cross-section. The result agrees with the consensus best values in the literature, and provides a derivation of k that does not rely on physical assumptions. & 2011 Elsevier Ltd. All rights reserved.
1. Introduction Timoshenko beam theory (TBT) provides shear deformation and rotatory inertia corrections to the classic Euler–Bernoulli theory [1]; it predicts the natural frequency of bending vibrations for long beams with remarkable accuracy if one employs the ‘‘best’’ value for the shear coefficient, k. Exact elastodynamic theory is available for beams of circular cross-section (Pochammer–Chree theory, see Love [2], article 202) and the thin (plane stress) rectangular section [3], and for these cases the best coefficients are k ¼ 6ð1 þ nÞ2 =ð7 þ 12n þ 4n2 Þ and k ¼ 5ð1þ nÞ=ð6 þ5nÞ, respectively, where n is Poisson’s ratio of the material. In turn, procedures have been developed for the general cross-section which lead to an expression for the best k in terms of the Saint-Venant flexure function, and which provide the above values when applied to these cross-sections. Stephen and Levinson [4,5] based their methods upon the static stress distribution for a beam subjected to gravity loading, rather than the tip loading assumed in the method proposed by Cowper [6]. More recently, Hutchinson [7] employed the Hellinger– Reissner variational principle to construct a beam theory of Timoshenko type, which incorporated an expression for the shear coefficient that was demonstrated to be equivalent to this best coefficient in the discussion and closure section of Ref. [7]. Hutchinson [8] provided further results for thin-walled beams. Despite these successes, all these works rely on ad hoc physical assumptions and are therefore sometimes queried. It would be of some advantage to be able to dispense with these assumptions, insightful though they are, and derive the shear coefficient for an arbitrary cross-section mathematically, and in the process reveal what approximations are in fact employed and therefore what the corrections are. In the present work, we consider standing waves in a beam of length L that is simply-supported (this restriction can be relaxed; see the end of Section 8.1); the governing elastodynamic equations are expanded as a power series in q ¼ mp=L, the integer m being the mode order; displacements, stress Corresponding author. Fax: + 852 2603 5252.
E-mail addresses: [email protected] (K.T. Chan), kfl[email protected] (K.F. Lai), [email protected] (N.G. Stephen), [email protected] (K. Young). 0022-460X/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jsv.2011.02.012
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
Nomenclature A e, E G i, j I J L m n q
cross-sectional area, constant dilatation, Young’s modulus shear modulus indices second moment of area integral related to Saint-Venant flexure function length of beam mode number normal, integer, order in q mp=L for a standing wave, wavenumber for a travelling wave
u, v, w x, y, z t, T U
g w k l
n r s, t o c
3489
displacement components Cartesian coordinates time, kinetic energy strain energy shear angle Saint-Venant flexure function shear coefficient Lame´ constant Poisson’s ratio density direct, shearing stress radian frequency cross-sectional rotation
components and frequencies are calculated for each power as necessary. The natural frequency is then expressed, for long thin beams (L large, q small) as a power series in q:
o2 ¼ A4 q4 þA6 q6 þ ,
(1)
in which symmetry only allows even powers of q. (All such series are meant to be asymptotic, not necessarily convergent; this is after all what is needed in applications with a fixed number of terms and q-0.) Euler–Bernoulli theory implies that the leading term is q4 and in fact gives the value of A4. The key in the present discussion is A6, and the strategy is to compute it in two different ways and compare the result. In Section 2, TBT is reviewed, and rendered into the form (1). The resultant A6 AT6 ðkÞ of course depends on k. Then we proceed with an alternate solution, by simply expanding the problem in powers of q, without relying on any physical assumptions or introducing any shear coefficient. The governing equations are set up in Section 3 and solved order-by-order in Section 4. The result is used to evaluate the strain energy U in Section 5 and the kinetic energy T in Section 6. The eigenvalue o2 is given by the Rayleigh quotient Q=U/T, which is evaluated in Section 7, giving a formula for A6 that does not contain k. Comparison with AT6 ðkÞ then yields k. The key result in (52) turns out to be identical with the canonical expression given by Stephen [4] and Stephen and Levinson [5] for an arbitrary cross-section, thus settling any possible controversy [9] that might remain. A discussion is given in Section 8 and a brief conclusion is given in Section 9. 2. Timoshenko beam theory We consider standing waves in a uniform, isotropic simply-supported beam of arbitrary cross-section and length L; the axial coordinate is z, and transverse vibration takes place in the xz-plane. Euler–Bernoulli theory considers just the transverse displacement u(z,t) and the curvature of the centre line. TBT expresses the centre-line slope in terms of the cross-sectional rotation cðz,tÞ and a centre-line shear angle gðz,tÞ ¼ cðz,tÞ þ quðz,tÞ=qz; the latter is related to the shear force by a shear coefficient k. Within this approximation, the coupled equations of free vibration may be written as q qu q2 u ¼ rA 2 , kAG cþ (2a) qz qz qt
kAG c þ
qu q2 c q2 c EIyy 2 ¼ rIyy 2 : qz qz qt
Elimination of c leads to a single fourth-order differential equation in both space and time 4 r2 Iyy q4 u q4 u q2 u E q u þ ¼ 0: EIyy 4 þ rA 2 rIyy 1 þ 2 2 kG qz qt kG qt 4 qz qt
(2b)
(3)
For a simply-supported beam the mode shape is sinusoidal in both space and time, so write upsinqzsinot,
(4)
where q ¼ mp=L is the wavenumber for mode m. The standing wave can be regarded as two superposed travelling waves, and strictly speaking the term ‘‘wavenumber’’ refers to the latter. Equivalently, one can also work with travelling waves, and use a complex notation, e.g., exp½iðqzotÞ, but for the present paper, the more physical notation of real variables will be used instead. One then has E r2 Iyy 4 q2 o2 þ EIyy q4 rAo2 rIyy 1 þ o ¼ 0: (5) kG kG
3490
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
Now the Euler–Bernoulli frequency oEB is defined by the first two terms in the above, that is,
o2EB ¼ ðEIyy =rAÞq4 , and it is convenient to embed this frequency into (5): divide throughout by EIyyq4 to give EI2yy 4 o 4 o 2 E Iyy 2 o 2 q q 1þ þ ¼ 0: 1 kG A oEB oEB oEB kA2 G
(6)
(7)
For long wavelengths, set all powers of q in (7) equal to zero, which leads to o ¼ oEB . As the wavelength becomes shorter, that is q becomes larger, so ðo=oEB Þ also becomes less than unity and one ignores the final term in (7), to give o 2 E Iyy 2 1 q ¼ 1þ 1 þ : (8) kG A oEB Using (6) for o2EB and employing the binomial expansion on the right-hand side of (8), one gets EI E Iyy 2 q þ Oðq8 Þ, o2 ¼ yy q4 1 1 þ rA kG A
(9)
which gives AT4 ¼
AT6 ¼
EIyy , rA
EIyy E Iyy 1þ , rA kG A
(10a)
(10b)
where the superscript T denotes that these come from TBT; AT4 merely expresses classical Euler–Bernoulli theory, while AT6 ðkÞ will allow k to be determined, if A6 can be found in an independent way—which will be the task of the rest of this paper. 3. Expansion of governing equations This section takes the fundamental equations of elastodynamics as applied to flexural vibrations, and simply performs an expansion in powers of q, with no other assumptions. The aim is to obtain (1) and in particular evaluate A6 without ever introducing k. This straightforward if apparently tedious task is made much easier by two important observations. First, the expression (1) may suggest the need to do a daunting sixth-order calculation in q. However, A4 and A6 can be evaluated if the eigenfunction is known up to third order (see Section 4.5). The ability to bypass the fourth-, fifth- and sixth-order eigenfunctions may appear somewhat fortuitous, but in fact exemplifies a general theorem [10]: if the eigenfunction is known with an error of O(qN) (here N= 4), and this is used in a Rayleigh quotient, then the eigenvalue can be evaluated with an error of O(q2N). Second, symmetry of the system under z-z, q-q implies that (1) involves only even powers of q, and more importantly, that the eigenfunctions are either even or odd in q, so that half the terms vanish. 3.1. Equations of motion Assume displacements of the form uðx,y,z,tÞ ¼ uðx,yÞsinqzsinot,
(11a)
vðx,y,z,tÞ ¼ vðx,yÞsinqzsinot,
(11b)
wðx,y,z,tÞ ¼ wðx,yÞcosqzsinot,
(11c)
where it is noticed that (u, v) and w are out of phase in z by a quarter cycle. The functions u, v and w are then expanded in powers of q. Symmetry implies that (u,v,w) must be odd in q, so u and v are even while w is odd in q: u ¼ u0 þu2 q2 þ ,
(12a)
v ¼ v0 þv2 q2 þ ,
(12b)
w ¼ w1 q þ w3 q3 þ :
(12c)
Terms beyond those shown are not necessary for our purpose. The dilatation is likewise expanded e
qu qv qw þ þ ¼ eðx,yÞsinqzsinot, qx qy qz
(13a)
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
e¼
qu0 qv0 qu2 qv2 þ þ þ w1 q2 þ , qx qy qx qy
3491
(13b)
where e is even in q. The Navier equations are three of the type ðl þGÞ 2
2
2
qe q2 u þGr2 u ¼ r 2 , qx qt
(14)
2
where r ¼ q =qx2 þ q =qy2 þ q =qz2 . Substituting the displacement (12) and (13) and the natural frequency according to (1), one then finds " # q2 q2 q 2 2 ðw ðv ðu þ u q þ Þ þ þ v q þ Þ q þ Þq ðl þGÞ 0 2 0 2 1 qx qxqy qx2 " # ! q2 q2 þ (15a) ðu0 þ u2 q2 þ Þðu0 þ u2 q2 þ Þq2 þ þ rðu0 þ u2 q2 þ ÞðA4 q4 þ A6 q6 þ Þ ¼ 0, þG qx2 qy2 "
" þG
# q2 q2 q 2 2 ðu0 þ u2 q þ Þ þ 2 ðv0 þ v2 q þ Þ ðw1 qþ Þq ðl þ GÞ qy qxqy qy # !
q2 q2 þ ðv0 þv2 q2 þ Þðv0 þv2 q2 þ Þq2 þ þ rðv0 þ v2 q2 þ ÞðA4 q4 þ A6 q6 þ Þ ¼ 0, qx2 qy2 q q ðu0 þ u2 q2 þ Þ þ ðv0 þ v2 q2 þ Þðw1 qþ Þq ðl þGÞq qx qy ! " # 2 2 q q 2 ðw þG þ q þ Þðw q þ Þq þ þ rðw1 qþ ÞðA4 q4 þ A6 q6 þ Þ ¼ 0, 1 1 qx2 qy2
(15b)
(15c)
which are to be solved order-by-order in q; even orders involve only (15a) and (15b), and odd orders only (15c). 3.2. Boundary conditions The boundary conditions of zero traction are
sx cosðx,nÞ þ txy cosðy,nÞ ¼ 0,
(16a)
txy cosðx,nÞ þ sy cosðy,nÞ ¼ 0,
(16b)
txz cosðx,nÞ þ tyz cosðy,nÞ ¼ 0,
(16c)
where si , i= x, y, z are the direct stresses and tij are the shear stresses; they also have a z- and t-dependence that can be factored out
si ðx,y,z,tÞ ¼ s i ðx,yÞsinqzsinot, i ¼ x,y,z,
(17a)
txy ðx,y,z,tÞ ¼ t xy ðx,yÞsinqzsinot,
(17b)
tiz ðx,y,z,tÞ ¼ t iz ðx,yÞcosqzsinot, i ¼ x,y:
(17c)
The direct and shear stress components are given in terms of the strains by
s x ¼ le þ 2Gðqu=qxÞ,
(18a)
s y ¼ le þ 2Gðqv=qyÞ,
(18b)
s z ¼ le2Gwq,
(18c)
t xy ¼ Gðqu=qy þ qv=qxÞ,
(18d)
t xz ¼ Gðuqþ qw=qxÞ,
(18e)
t yz ¼ Gðvqþ qw=qyÞ:
(18f)
s x ¼ le þ 2G½ðqu0 =qxÞ þðqu2 =qxÞq2 þ ,
(19a)
s y ¼ le þ 2G½ðqv0 =qyÞ þ ðqv2 =qyÞq2 þ ,
(19b)
In more detail,
3492
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
s z ¼ le2G½w1 q2 þ w3 q4 þ ,
(19c)
t xy ¼ G½ðqu0 =qy þ qv0 =qxÞ þ ðqu2 =qyþ qv2 =qxÞq2 þ ,
(19d)
t xz ¼ G½ðu0 þ qw1 =qxÞqþ ,
(19e)
t yz ¼ G½ðv0 þ qw1 =qyÞq þ :
(19f)
We will also find it convenient to write these components as
s i ¼ sð0Þ þ sð2Þ q2 þ , i ¼ x,y,z, i i
(20a)
ð0Þ 2 t xy ¼ txy þ tð2Þ xy q þ ,
(20b)
ð3Þ 3 t iz ¼ tð1Þ i ¼ x,y, iz q þ tiz q þ ,
(20c)
in which all terms with the wrong symmetry in q have been dropped. (If these were kept at this point, we would simply find, upon calculation order-by-order, that they in fact vanish.) 4. Order-by-order solution In this section, we solve the equations of motion (15) subject to the boundary conditions (16), order-by-order. 4.1. Zeroth order Suppose, as a matter of convention, that the lowest order displacement is in the x-direction, with the normalization 2 set to unity. (Otherwise, all amplitudes will carry a factor u0 and all energies a factor u0, which will in the end cancel in Q= U/T.) Thus u0 ¼ 1,
v0 ¼ 0
(21)
with w0 = 0 already assumed in (12). It is obvious that the equations of motion (15) are satisfied, and since all stress components are zero to this order, the boundary conditions (16) are obviously satisfied as well. 4.2. First order The solution is obviously w1 ¼ x
(22)
with u1 = v1 = 0 already assumed in (12). The Navier equations (15) are satisfied, and
tð1Þ xz ¼ Gðu0 þ qw1 =qxÞ ¼ 0
(23)
by (21) and (22). All other stress components are obviously zero as well. Again, the boundary conditions need not be considered. For later reference, it is important to note that the stresses, therefore, start at O(q2), and hence the strain energy at O(q4). 4.3. Second order The solution is u2 ¼ nðy2 x2 Þ=2,
v2 ¼ nxy
(24)
with w2 = 0 already assumed in (12). The only non-zero stress is
sð2Þ z ¼ Ex:
(25)
This is equivalent to the Euler–Bernoulli stress distribution. 4.4. Third order For the third order, we only need to determine w3, since symmetry dictates u3 =v3 = 0. Only (15c) needs to be considered, and this reduces to
r2 w3 ¼ 2x,
(26)
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497 2
3493 2
where henceforth and without danger of confusion r stands for the two-dimensional Laplacian: r ¼ q2 =qx2 þq2 =qy2 . Guided by Love [2], set w3 ¼ ðw þ xy2 Þ,
(27)
r2 w ¼ 0:
(28)
so that
To determine the boundary condition on w, we first note that the two non-zero third-order stresses are 2 2 tð3Þ xz ¼ Gðu2 þ qw3 =qxÞ ¼ G½qw=qx þ nx =2þ ð1n=2Þy ,
(29a)
tð3Þ yz ¼ Gðv2 þ qw3 =qyÞ ¼ G½qw=qy þ ð2 þ nÞxy:
(29b)
We only need to consider (16c), which reduces to qw n 2 n qw þ x þ 1 y2 cosðx,nÞ þ þ ð2þ nÞy2 cosðy,nÞ ¼ 0: qx 2 2 qy
(30)
But qw qw qw dx qw dy dw þ ¼ cosðx,nÞ þ cosðy,nÞ ¼ , qx qy qx dn qy dn dn
(31)
hn dw n i ¼ x2 þ 1 y2 cosðx,nÞð2 þ nÞxycosðy,nÞ: dn 2 2
(32)
so that (30) becomes
Thus w is determined by the differential equation (28) together with the Neumann boundary condition (32), and is seen to be nothing other than the Saint-Venant flexure function; see Ref. [2, Chapter XV]. 4.5. Higher orders As far as the eigenvalue to O(q6) is concerned, the fourth and higher order eigenfunctions will cancel when put into the Rayleigh quotient. This ‘‘miraculous’’ cancellation is best understood in a general context [10]. But the upshot for the present purpose is that these higher order eigenfunctions need not be evaluated. 5. Strain energy In terms of the stress components, the strain energy of the beam is given by Ref. [11], article 90: ZZ L 1 n 1 2 ðs 2x þ s 2y þ s 2z Þ ðs x s y þ s y s z þ s z s x Þ þ ðt xy þ t 2yz þ t 2zx Þ dx dy, U¼ 2 2E 2G E
(33)
in which a trivial integration over z has been carried out. Denote U ¼ U0 þ U2 q2 þ :
(34)
5.1. Zeroth order ð0Þ2 This will involve terms such as txy through to szð0Þ2 , as well as terms such as sxð0Þ sð0Þ y , but since all zeroth-order stress components are zero, one immediately has U0 = 0.
5.2. First order ð1Þ ð0Þ ð1Þ ð0Þ ð1Þ ð1Þ ð0Þ This will involve terms such as 2tð0Þ xy txy through to 2sz sz , as well as terms such as sx sy þ sx sy , but since all zeroth-order and first-order stress components are zero, one immediately has U1 =0.
5.3. Second order ð0Þ ð2Þ ð1Þ2 ð0Þ ð2Þ This will involve terms such as tð1Þ2 þ2szð0Þ sð2Þ xy þ2txy txy through to sz z , as well as terms such as sx sy þ ð1Þ ð2Þ ð0Þ s sy þ sx sy ; again all zeroth-order and first-order stress components are zero, and one immediately has U2 =0. ð1Þ x
5.4. Third order ð3Þ ð1Þ ð2Þ ð0Þ ð3Þ ð1Þ ð2Þ ð0Þ ð3Þ This will involve terms such as 2tð0Þ xy txy þ 2txy txy through to 2sz sz þ 2sz sz , as well as terms such as sx sy þ ð2Þ ð2Þ ð1Þ ð3Þ ð0Þ sð1Þ s þ s s þ s s . Again one immediately has U = 0. 3 x y x y x y
3494
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
5.5. Fourth order ð0Þ ð4Þ ð1Þ ð3Þ ð2Þ2 ð4Þ ð1Þ ð3Þ This will involve terms such as tð2Þ2 þ 2sð0Þ xy þ 2txy txy þ 2txy txy through to sz z sz þ 2sz sz , as well as terms such as ð4Þ ð1Þ ð3Þ ð2Þ ð2Þ ð3Þ ð1Þ ð4Þ ð0Þ sð0Þ x sy þ sx sy þ sx sy þ sx sy þ sx sy . The only non-zero contributor is
sð2Þ z ¼ Ex,
(35)
U4 ¼ EIyy L=4:
(36)
and one finds
5.6. Fifth order ð0Þ ð5Þ ð4Þ ð2Þ ð3Þ ð0Þ ð5Þ ð1Þ ð4Þ ð2Þ ð3Þ txy þ2tð1Þ This involves terms such as 2txy xy txy þ 2txy txy through to 2sz sz þ2sz sz þ 2sz sz , as well as terms such ð5Þ ð1Þ ð4Þ ð2Þ ð3Þ ð3Þ ð2Þ ð4Þ ð1Þ ð5Þ ð0Þ as sð0Þ x sy þ sx sy þ sx sy þ sx sy þ sx sy þ sx sy . The only possible contributors come from the second- and thirdð3Þ ð3Þ order stresses, but since there is no product involving sð2Þ z and txz or tyz , one immediately has U5 = 0. This result (and the same for U1 and U3) is anticipated since the strain energy must be even in q.
5.7. Sixth order ð3Þ2 At this level, there are a variety of terms which do contribute; these are tð3Þ2 xz and tyz , which are straightforward, and ð2Þ ð4Þ ð2Þ ð4Þ ð2Þ ð4Þ also 2sð4Þ z sz , sx sz and sy sz . This suggests that one must determine the direct stress components si , i= x, y, z, but in fact one only needs a knowledge of w3. In terms of the displacement components, one has ð4Þ sð4Þ þ2G x ¼ le
qu4 , qx
(37a)
ð4Þ sð4Þ þ2G y ¼ le
qv4 , qy
(37b)
ð4Þ sð4Þ z ¼ le 2Gw3 ,
(37c)
where eð4Þ ¼
qu4 qv4 þ w3 : qx qy
(38)
The relevant expression in the integrand in U is 1 n ð4Þ ð2Þ ð2Þ ð4Þ ð2Þ ð2sð4Þ z sz Þ ðsy sz þ sy sz Þ, 2E E
(39)
sð2Þ z w3 ,
(40)
which reduces to
where w3 comes from the last term in (19c) and all reference to eigenfunctions beyond third order has disappeared. One then finds EL 1 EL J1 þ J, (41) J2 U6 ¼ 2 4ð1 þ nÞ 2 in terms of two integrals with dimensions of (length)6 defined in terms of the Saint-Venant flexure function: ZZ J1 ¼ xðw þ xy2 Þ dx dy, 2 ) ZZ ( qw n 2 n 2 2 qw þ x þ 1 y þð2 þ nÞxy J2 ¼ þ dx dy: qx 2 2 qy
(42a)
(42b)
6. Kinetic energy The Rayleigh quotient Q (see next section) is essentially the ratio of the strain energy U to the kinetic energy T. (This term is used as a shorthand. The kinetic energy is actually o2 T.) Since U starts with q4 and we want Q only to q6, it suffices to calculate T to just q2. Now we have ZZ rL ðu 2 þv 2 þ w 2 Þ dx dy, T¼ (43) 4
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
3495
where again the trivial integration over z has been carried out. The integrand is ðu0 þ u2 q2 þ Þ2 þ ðv0 þv2 q2 þ Þ2 þðw1 qþ w3 q3 Þ2 : The integral of
2 u0
0
provides the q term, while
2 w1 +2u0u2
(44)
2
provides the q term, and one finds
T ¼ T0 þT2 q2 þ ,
(45)
where
T2 ¼
rL 4
rAL
,
(46a)
½Iyy þ nðIxx Iyy Þ:
(46b)
T0 ¼
4
7. Rayleigh quotient It is well known that the eigenvalue o2 is given by the Rayleigh quotient [12]
o2 ¼ Q
U U4 q4 þ U6 q6 þ ¼ , T T0 þT2 q2 þ
(47)
where U and T are evaluated for the corresponding eigenfunction; from this one finds the coefficients defined by (1) as U4 EIyy , ¼ T0 rA
(48a)
U6 U4 T2 E 2 ¼ f2AJIyy ½Iyy þ nðIxx Iyy Þg, T0 rA2 T0
(48b)
A4 ¼
A6 ¼ so that compared with (10b), we find
k¼
2 2ð1þ nÞIyy
2AJ 1 þ AJ 2 =½2ð1 þ nÞ þ nIyy ðIyy Ixx Þ
:
(49)
This expression is similar but not identical to that given by Hutchinson [7]. But using an identity presented in the discussion of Ref. [7], we have J2 ¼ nJ3 2ð1 þ nÞJ1 ,
(50)
where J3 ¼
ZZ 2 2 x y qw n 2 n qw þ x þ 1 y2 þ xy þ ð2 þ nÞxy dx dy: 2 qx 2 2 qy
(51)
Thus we can finally render (49) into
k¼
2 4ð1þ nÞ2 Iyy : 2ð1 þ nÞAJ 1 þ nAJ 3 þ 2nð1 þ nÞIyy ðIyy Ixx Þ
(52)
In this form, k agrees exactly with the expression for the shear coefficient presented in Refs. [4,5], thus proving the latter without having to resort to the physical assumption of TBT. 8. Discussion 8.1. General remarks Our derivation relies on a single approximation, namely that the wavelength is long, i.e., q is small, so that a power series in q makes sense. There is no need to guess, on physical grounds, what degrees of freedom must be kept. In fact, the new variable g in TBT emerges automatically, in the following way. First, the centre line of the beam tilts by an angle (dropping common time-dependent factors from (11)) qu ¼ qucosqz ¼ qðu0 þ u2 q2 þ Þcosqz: qz On the other hand, the cross-sectional plane tilts by an angle qw qw qw1 qw3 3 ¼ cosqz ¼ qþ q þ cosqz: c¼ qx qx qx qx
(53)
(54)
3496
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
The centre line and the cross-sectional plane will, therefore, deviate from orthogonality by an angle qu qw1 qw3 3 ¼ u0 þ q þ u2 þ q þ cosqz: g ¼ cþ qz qx qx
(55)
Our solution, with no assumptions, shows that the first bracket is zero—recovering the key insight and assumption g ¼ 0 in Euler–Bernoulli theory. The next bracket is not zero, so g must be kept for shorter wavelengths, and moreover its effect is captured once w3 is evaluated, without the need to introduce any parameters. In principle, the expansion in q can be continued and is formally exact, although the evaluation of u4, v4, w5, yis impossible in practice—with two important exceptions to be discussed below. But some general features of the expansion can be noted. First, purely on dimensional grounds, each successive term in (1) has an extra power of (qa)2, where a is a typical transverse dimension of the beam, i.e., A2n ðE=ra2 ÞðqaÞ2n . Second, we do not expect TBT to reproduce the next order exactly, i.e., AT8 aA8 , since for shorter wavelengths, the vibration must be described by more than two variables uðz,tÞ, cðz,tÞ at each z — as can be demonstrated in simple cases (see below). These remarks, taken together, imply that TBT, though highly accurate when qa 5 1, cannot be expected to work when qa 1, since the A8q8yterms become important and cannot be reproduced correctly. The discussion in this paper refers to a simply-supported beam, i.e., hinged–hinged end conditions, but this restriction is unnecessary, since other conditions, e.g., guided–guided or guided–hinged, can be regarded as portions of a multi-span hinged–hinged beam [13]. It should also be mentioned that our ‘‘best’’ choice of k is obtained by matching the q6 term in the dispersion relation of the lowest branch (the one without nodes in the cross-sectional plane), and is of course the optimal one for using TBT to describe oscillations of this type—which are the ones most commonly encountered in engineering practice. If one were interested in other types of oscillations, for example the higher branches, other choices may be more appropriate. 8.2. Solvable examples and possible generalizations In two cases, the expansion in q can actually be carried out to very high orders. A brief discussion is given here, principally to illustrate the qualitative remarks above in a precise setting. First consider flexural vibrations in a hypothetical world of two dimensions, say xz. The partial differential equations in Sections 3 and 4 become ordinary differential equations in x. Moreover, on dimensional grounds, the n th-order eigenfunction must go as ðqaÞn ðx=aÞk , with k r n to ensure regularity when a-0; thus it must be a polynomial in x of maximum order n. The differential equation can be cast into algebraic recursion relations for the polynomial coefficients. With these simplifications, the solution can be carried out to many orders. Next consider longitudinal vibrations in a circular cylinder. Using cylindrical coordinates ðr, y,zÞ and factoring out the trivial y-dependence, one again obtains an ordinary differential equation (though of a slightly more complicated form) in r and for the same dimensional reason, the solution is again a polynomial in r, which likewise can be found to very high orders. We have solved both of these cases in powers of q to 20 orders. The coefficients in the polynomials turn out to be rational functions of n, involving powers up to n3n3 ; though complicated, these can be obtained using algebraic software packages. If one seeks only numerical values for a specific value of n, the computation is much simpler. The results, which can also be checked against exact solutions available in these cases [14] confirm the qualitative features discussed in the last subsection. Incidentally, an extension of the two-dimensional problem is a thin plates of thickness h, which may be treated in a manner similar to that given in the present paper, to provide a derivation of the ‘‘best’’ value of k for the Mindlin theory of thin plates [15]. That study will be given elsewhere. 8.3. Using the Rayleigh quotient The present paper makes use of the Rayleigh quotient, which has the advantage [10] of giving the sixth-order eigenvalue A6 from the third-order eigenfunction w3. Two other nice features should be mentioned as well. First, using any approximate eigenfunction to evaluate U and T, and hence Q, guarantees positivity. In contrast, the approximation o2 ¼ A4 q4 þA6 q6 goes negative for qa = O(1). Second, again on physical grounds, one expects o2 =q2 - constant as q-1 (finite phase velocity). This property is also nicely guaranteed for the Rayleigh quotient, since U involves two extra powers of q compared to T. Because of these properties, the Rayleigh quotient is often accurate over a wider range of qa. For longitudinal vibrations in a circular cylinder, this method leads to 5 percent accuracy up to qa 2:5, namely a cylinder with diameter larger than its length—almost a ‘‘disk’’ rather than a ‘‘rod’’. 9. Conclusion In conclusion, the simple and straightforward strategy of expanding in powers of q provides an alternative method to evaluate the shear coefficient k that is systematic and unambiguous. The expansion, formally exact when carried out to all
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
3497
orders, also provides a wider perspective to view TBT: (a) with the ‘‘best’’ value for k, it gives the next OððqaÞ2 Þ correction to the classical Euler–Bernoulli theory, but (b) it is itself unlikely to be accurate when qa 1. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
S.P. Timoshenko, On the transverse vibrations of bars of uniform cross-section, Philosophical Magazine 43 (1922) 125–131. A.E.H. Love, A Treatise on the Mathematical Theory of Elasticity, fourth ed., Dover, New York, 1944. N.G. Stephen, On ‘‘A check on the accuracy of Timoshenko’s beam theory’’, Journal of Sound and Vibration 257 (2002) 809–812. N.G. Stephen, Timoshenko’s shear coefficient from a beam subjected to gravity loading, ASME Journal of Applied Mechanics 47 (1980) 121–127. N.G. Stephen, M. Levinson, A second order beam theory, Journal of Sound and Vibration 67 (1979) 293–305. G.R. Cowper, The shear coefficient in Timoshenko’s beam theory, ASME Journal of Applied Mechanics 33 (1966) 335–340. J.R. Hutchinson, Shear coefficients for Timoshenko beam theory, ASME Journal of Applied Mechanics 68 (2001) 87–92 See also Discussion and Closure, ASME Journal of Applied Mechanics 68 (2001) 959–961. J.R. Hutchinson, Shear coefficients for thin-walled Timoshenko beams, Third International Symposium on the Vibrations of Continuous Systems, Jackson Hole, Wyoming, USA, July 2001. J.D. Renton, A check on the accuracy of Timoshenko’s beam theory, Journal of Sound and Vibration 245 (2001) 559–561. K.T. Chan, N.G. Stephen, K. Young, Perturbation theory and the Rayleigh quotient. Journal of Sound and Vibration, 330 (2011) 2073–2078 (Ref: Ms. No. JSVD-10-00949R1). S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, third ed., McGraw-Hill, 1970. S.S. Rao, Mechanical Vibrations, fourth ed., Pearson Prentice Hall, Upper Saddle River, NJ, USA, 2004. N.G. Stephen, The second spectrum of Timoshenko beam theory—further assessment, Journal of Sound and Vibration 292 (2006) 372–389. A.E. Armenakis, D.C. Gazis, G. Herrmann, Free Vibration of Circular Cylindrical Shells, Pergamon, Oxford, 1969. R.D. Mindlin, Influence of rotatory inertia and shear on flexural motions of isotropic elastic plates, ASME Journal of Applied Mechanics 18 (1951) 31–38.
Contents lists available at ScienceDirect
Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi
A new method to determine the shear coefficient of Timoshenko beam theory K.T. Chan a, K.F. Lai b, N.G. Stephen c, K. Young b, a b c
Department of Mechanical Engineering, The Hong Kong Polytechnic University, Hong Kong, PR China Department of Physics, The Chinese University of Hong Kong, Hong Kong, PR China School of Engineering Sciences, Mechanical Engineering, University of Southampton, Highfield, Southampton SO17 1BJ, UK
a r t i c l e in f o
abstract
Article history: Received 17 November 2010 Received in revised form 8 February 2011 Accepted 8 February 2011 Handling Editor: S. Ilanko Available online 5 March 2011
The frequencies o of flexural vibrations in a uniform beam of arbitrary cross-section and length L are analysed by expanding the exact elastodynamics equations in powers of the wavenumber q ¼ mp=L, where m is the mode number: o2 ¼ A4 q4 þ A6 q6 þ . The coefficients A4 and A6 are obtained without further assumptions; the former captures Euler–Bernoulli theory while the latter, when compared with Timoshenko beam theory rendered into the same form, unambiguously yields the shear coefficient k for any cross-section. The result agrees with the consensus best values in the literature, and provides a derivation of k that does not rely on physical assumptions. & 2011 Elsevier Ltd. All rights reserved.
1. Introduction Timoshenko beam theory (TBT) provides shear deformation and rotatory inertia corrections to the classic Euler–Bernoulli theory [1]; it predicts the natural frequency of bending vibrations for long beams with remarkable accuracy if one employs the ‘‘best’’ value for the shear coefficient, k. Exact elastodynamic theory is available for beams of circular cross-section (Pochammer–Chree theory, see Love [2], article 202) and the thin (plane stress) rectangular section [3], and for these cases the best coefficients are k ¼ 6ð1 þ nÞ2 =ð7 þ 12n þ 4n2 Þ and k ¼ 5ð1þ nÞ=ð6 þ5nÞ, respectively, where n is Poisson’s ratio of the material. In turn, procedures have been developed for the general cross-section which lead to an expression for the best k in terms of the Saint-Venant flexure function, and which provide the above values when applied to these cross-sections. Stephen and Levinson [4,5] based their methods upon the static stress distribution for a beam subjected to gravity loading, rather than the tip loading assumed in the method proposed by Cowper [6]. More recently, Hutchinson [7] employed the Hellinger– Reissner variational principle to construct a beam theory of Timoshenko type, which incorporated an expression for the shear coefficient that was demonstrated to be equivalent to this best coefficient in the discussion and closure section of Ref. [7]. Hutchinson [8] provided further results for thin-walled beams. Despite these successes, all these works rely on ad hoc physical assumptions and are therefore sometimes queried. It would be of some advantage to be able to dispense with these assumptions, insightful though they are, and derive the shear coefficient for an arbitrary cross-section mathematically, and in the process reveal what approximations are in fact employed and therefore what the corrections are. In the present work, we consider standing waves in a beam of length L that is simply-supported (this restriction can be relaxed; see the end of Section 8.1); the governing elastodynamic equations are expanded as a power series in q ¼ mp=L, the integer m being the mode order; displacements, stress Corresponding author. Fax: + 852 2603 5252.
E-mail addresses: [email protected] (K.T. Chan), kfl[email protected] (K.F. Lai), [email protected] (N.G. Stephen), [email protected] (K. Young). 0022-460X/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jsv.2011.02.012
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
Nomenclature A e, E G i, j I J L m n q
cross-sectional area, constant dilatation, Young’s modulus shear modulus indices second moment of area integral related to Saint-Venant flexure function length of beam mode number normal, integer, order in q mp=L for a standing wave, wavenumber for a travelling wave
u, v, w x, y, z t, T U
g w k l
n r s, t o c
3489
displacement components Cartesian coordinates time, kinetic energy strain energy shear angle Saint-Venant flexure function shear coefficient Lame´ constant Poisson’s ratio density direct, shearing stress radian frequency cross-sectional rotation
components and frequencies are calculated for each power as necessary. The natural frequency is then expressed, for long thin beams (L large, q small) as a power series in q:
o2 ¼ A4 q4 þA6 q6 þ ,
(1)
in which symmetry only allows even powers of q. (All such series are meant to be asymptotic, not necessarily convergent; this is after all what is needed in applications with a fixed number of terms and q-0.) Euler–Bernoulli theory implies that the leading term is q4 and in fact gives the value of A4. The key in the present discussion is A6, and the strategy is to compute it in two different ways and compare the result. In Section 2, TBT is reviewed, and rendered into the form (1). The resultant A6 AT6 ðkÞ of course depends on k. Then we proceed with an alternate solution, by simply expanding the problem in powers of q, without relying on any physical assumptions or introducing any shear coefficient. The governing equations are set up in Section 3 and solved order-by-order in Section 4. The result is used to evaluate the strain energy U in Section 5 and the kinetic energy T in Section 6. The eigenvalue o2 is given by the Rayleigh quotient Q=U/T, which is evaluated in Section 7, giving a formula for A6 that does not contain k. Comparison with AT6 ðkÞ then yields k. The key result in (52) turns out to be identical with the canonical expression given by Stephen [4] and Stephen and Levinson [5] for an arbitrary cross-section, thus settling any possible controversy [9] that might remain. A discussion is given in Section 8 and a brief conclusion is given in Section 9. 2. Timoshenko beam theory We consider standing waves in a uniform, isotropic simply-supported beam of arbitrary cross-section and length L; the axial coordinate is z, and transverse vibration takes place in the xz-plane. Euler–Bernoulli theory considers just the transverse displacement u(z,t) and the curvature of the centre line. TBT expresses the centre-line slope in terms of the cross-sectional rotation cðz,tÞ and a centre-line shear angle gðz,tÞ ¼ cðz,tÞ þ quðz,tÞ=qz; the latter is related to the shear force by a shear coefficient k. Within this approximation, the coupled equations of free vibration may be written as q qu q2 u ¼ rA 2 , kAG cþ (2a) qz qz qt
kAG c þ
qu q2 c q2 c EIyy 2 ¼ rIyy 2 : qz qz qt
Elimination of c leads to a single fourth-order differential equation in both space and time 4 r2 Iyy q4 u q4 u q2 u E q u þ ¼ 0: EIyy 4 þ rA 2 rIyy 1 þ 2 2 kG qz qt kG qt 4 qz qt
(2b)
(3)
For a simply-supported beam the mode shape is sinusoidal in both space and time, so write upsinqzsinot,
(4)
where q ¼ mp=L is the wavenumber for mode m. The standing wave can be regarded as two superposed travelling waves, and strictly speaking the term ‘‘wavenumber’’ refers to the latter. Equivalently, one can also work with travelling waves, and use a complex notation, e.g., exp½iðqzotÞ, but for the present paper, the more physical notation of real variables will be used instead. One then has E r2 Iyy 4 q2 o2 þ EIyy q4 rAo2 rIyy 1 þ o ¼ 0: (5) kG kG
3490
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
Now the Euler–Bernoulli frequency oEB is defined by the first two terms in the above, that is,
o2EB ¼ ðEIyy =rAÞq4 , and it is convenient to embed this frequency into (5): divide throughout by EIyyq4 to give EI2yy 4 o 4 o 2 E Iyy 2 o 2 q q 1þ þ ¼ 0: 1 kG A oEB oEB oEB kA2 G
(6)
(7)
For long wavelengths, set all powers of q in (7) equal to zero, which leads to o ¼ oEB . As the wavelength becomes shorter, that is q becomes larger, so ðo=oEB Þ also becomes less than unity and one ignores the final term in (7), to give o 2 E Iyy 2 1 q ¼ 1þ 1 þ : (8) kG A oEB Using (6) for o2EB and employing the binomial expansion on the right-hand side of (8), one gets EI E Iyy 2 q þ Oðq8 Þ, o2 ¼ yy q4 1 1 þ rA kG A
(9)
which gives AT4 ¼
AT6 ¼
EIyy , rA
EIyy E Iyy 1þ , rA kG A
(10a)
(10b)
where the superscript T denotes that these come from TBT; AT4 merely expresses classical Euler–Bernoulli theory, while AT6 ðkÞ will allow k to be determined, if A6 can be found in an independent way—which will be the task of the rest of this paper. 3. Expansion of governing equations This section takes the fundamental equations of elastodynamics as applied to flexural vibrations, and simply performs an expansion in powers of q, with no other assumptions. The aim is to obtain (1) and in particular evaluate A6 without ever introducing k. This straightforward if apparently tedious task is made much easier by two important observations. First, the expression (1) may suggest the need to do a daunting sixth-order calculation in q. However, A4 and A6 can be evaluated if the eigenfunction is known up to third order (see Section 4.5). The ability to bypass the fourth-, fifth- and sixth-order eigenfunctions may appear somewhat fortuitous, but in fact exemplifies a general theorem [10]: if the eigenfunction is known with an error of O(qN) (here N= 4), and this is used in a Rayleigh quotient, then the eigenvalue can be evaluated with an error of O(q2N). Second, symmetry of the system under z-z, q-q implies that (1) involves only even powers of q, and more importantly, that the eigenfunctions are either even or odd in q, so that half the terms vanish. 3.1. Equations of motion Assume displacements of the form uðx,y,z,tÞ ¼ uðx,yÞsinqzsinot,
(11a)
vðx,y,z,tÞ ¼ vðx,yÞsinqzsinot,
(11b)
wðx,y,z,tÞ ¼ wðx,yÞcosqzsinot,
(11c)
where it is noticed that (u, v) and w are out of phase in z by a quarter cycle. The functions u, v and w are then expanded in powers of q. Symmetry implies that (u,v,w) must be odd in q, so u and v are even while w is odd in q: u ¼ u0 þu2 q2 þ ,
(12a)
v ¼ v0 þv2 q2 þ ,
(12b)
w ¼ w1 q þ w3 q3 þ :
(12c)
Terms beyond those shown are not necessary for our purpose. The dilatation is likewise expanded e
qu qv qw þ þ ¼ eðx,yÞsinqzsinot, qx qy qz
(13a)
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
e¼
qu0 qv0 qu2 qv2 þ þ þ w1 q2 þ , qx qy qx qy
3491
(13b)
where e is even in q. The Navier equations are three of the type ðl þGÞ 2
2
2
qe q2 u þGr2 u ¼ r 2 , qx qt
(14)
2
where r ¼ q =qx2 þ q =qy2 þ q =qz2 . Substituting the displacement (12) and (13) and the natural frequency according to (1), one then finds " # q2 q2 q 2 2 ðw ðv ðu þ u q þ Þ þ þ v q þ Þ q þ Þq ðl þGÞ 0 2 0 2 1 qx qxqy qx2 " # ! q2 q2 þ (15a) ðu0 þ u2 q2 þ Þðu0 þ u2 q2 þ Þq2 þ þ rðu0 þ u2 q2 þ ÞðA4 q4 þ A6 q6 þ Þ ¼ 0, þG qx2 qy2 "
" þG
# q2 q2 q 2 2 ðu0 þ u2 q þ Þ þ 2 ðv0 þ v2 q þ Þ ðw1 qþ Þq ðl þ GÞ qy qxqy qy # !
q2 q2 þ ðv0 þv2 q2 þ Þðv0 þv2 q2 þ Þq2 þ þ rðv0 þ v2 q2 þ ÞðA4 q4 þ A6 q6 þ Þ ¼ 0, qx2 qy2 q q ðu0 þ u2 q2 þ Þ þ ðv0 þ v2 q2 þ Þðw1 qþ Þq ðl þGÞq qx qy ! " # 2 2 q q 2 ðw þG þ q þ Þðw q þ Þq þ þ rðw1 qþ ÞðA4 q4 þ A6 q6 þ Þ ¼ 0, 1 1 qx2 qy2
(15b)
(15c)
which are to be solved order-by-order in q; even orders involve only (15a) and (15b), and odd orders only (15c). 3.2. Boundary conditions The boundary conditions of zero traction are
sx cosðx,nÞ þ txy cosðy,nÞ ¼ 0,
(16a)
txy cosðx,nÞ þ sy cosðy,nÞ ¼ 0,
(16b)
txz cosðx,nÞ þ tyz cosðy,nÞ ¼ 0,
(16c)
where si , i= x, y, z are the direct stresses and tij are the shear stresses; they also have a z- and t-dependence that can be factored out
si ðx,y,z,tÞ ¼ s i ðx,yÞsinqzsinot, i ¼ x,y,z,
(17a)
txy ðx,y,z,tÞ ¼ t xy ðx,yÞsinqzsinot,
(17b)
tiz ðx,y,z,tÞ ¼ t iz ðx,yÞcosqzsinot, i ¼ x,y:
(17c)
The direct and shear stress components are given in terms of the strains by
s x ¼ le þ 2Gðqu=qxÞ,
(18a)
s y ¼ le þ 2Gðqv=qyÞ,
(18b)
s z ¼ le2Gwq,
(18c)
t xy ¼ Gðqu=qy þ qv=qxÞ,
(18d)
t xz ¼ Gðuqþ qw=qxÞ,
(18e)
t yz ¼ Gðvqþ qw=qyÞ:
(18f)
s x ¼ le þ 2G½ðqu0 =qxÞ þðqu2 =qxÞq2 þ ,
(19a)
s y ¼ le þ 2G½ðqv0 =qyÞ þ ðqv2 =qyÞq2 þ ,
(19b)
In more detail,
3492
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
s z ¼ le2G½w1 q2 þ w3 q4 þ ,
(19c)
t xy ¼ G½ðqu0 =qy þ qv0 =qxÞ þ ðqu2 =qyþ qv2 =qxÞq2 þ ,
(19d)
t xz ¼ G½ðu0 þ qw1 =qxÞqþ ,
(19e)
t yz ¼ G½ðv0 þ qw1 =qyÞq þ :
(19f)
We will also find it convenient to write these components as
s i ¼ sð0Þ þ sð2Þ q2 þ , i ¼ x,y,z, i i
(20a)
ð0Þ 2 t xy ¼ txy þ tð2Þ xy q þ ,
(20b)
ð3Þ 3 t iz ¼ tð1Þ i ¼ x,y, iz q þ tiz q þ ,
(20c)
in which all terms with the wrong symmetry in q have been dropped. (If these were kept at this point, we would simply find, upon calculation order-by-order, that they in fact vanish.) 4. Order-by-order solution In this section, we solve the equations of motion (15) subject to the boundary conditions (16), order-by-order. 4.1. Zeroth order Suppose, as a matter of convention, that the lowest order displacement is in the x-direction, with the normalization 2 set to unity. (Otherwise, all amplitudes will carry a factor u0 and all energies a factor u0, which will in the end cancel in Q= U/T.) Thus u0 ¼ 1,
v0 ¼ 0
(21)
with w0 = 0 already assumed in (12). It is obvious that the equations of motion (15) are satisfied, and since all stress components are zero to this order, the boundary conditions (16) are obviously satisfied as well. 4.2. First order The solution is obviously w1 ¼ x
(22)
with u1 = v1 = 0 already assumed in (12). The Navier equations (15) are satisfied, and
tð1Þ xz ¼ Gðu0 þ qw1 =qxÞ ¼ 0
(23)
by (21) and (22). All other stress components are obviously zero as well. Again, the boundary conditions need not be considered. For later reference, it is important to note that the stresses, therefore, start at O(q2), and hence the strain energy at O(q4). 4.3. Second order The solution is u2 ¼ nðy2 x2 Þ=2,
v2 ¼ nxy
(24)
with w2 = 0 already assumed in (12). The only non-zero stress is
sð2Þ z ¼ Ex:
(25)
This is equivalent to the Euler–Bernoulli stress distribution. 4.4. Third order For the third order, we only need to determine w3, since symmetry dictates u3 =v3 = 0. Only (15c) needs to be considered, and this reduces to
r2 w3 ¼ 2x,
(26)
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497 2
3493 2
where henceforth and without danger of confusion r stands for the two-dimensional Laplacian: r ¼ q2 =qx2 þq2 =qy2 . Guided by Love [2], set w3 ¼ ðw þ xy2 Þ,
(27)
r2 w ¼ 0:
(28)
so that
To determine the boundary condition on w, we first note that the two non-zero third-order stresses are 2 2 tð3Þ xz ¼ Gðu2 þ qw3 =qxÞ ¼ G½qw=qx þ nx =2þ ð1n=2Þy ,
(29a)
tð3Þ yz ¼ Gðv2 þ qw3 =qyÞ ¼ G½qw=qy þ ð2 þ nÞxy:
(29b)
We only need to consider (16c), which reduces to qw n 2 n qw þ x þ 1 y2 cosðx,nÞ þ þ ð2þ nÞy2 cosðy,nÞ ¼ 0: qx 2 2 qy
(30)
But qw qw qw dx qw dy dw þ ¼ cosðx,nÞ þ cosðy,nÞ ¼ , qx qy qx dn qy dn dn
(31)
hn dw n i ¼ x2 þ 1 y2 cosðx,nÞð2 þ nÞxycosðy,nÞ: dn 2 2
(32)
so that (30) becomes
Thus w is determined by the differential equation (28) together with the Neumann boundary condition (32), and is seen to be nothing other than the Saint-Venant flexure function; see Ref. [2, Chapter XV]. 4.5. Higher orders As far as the eigenvalue to O(q6) is concerned, the fourth and higher order eigenfunctions will cancel when put into the Rayleigh quotient. This ‘‘miraculous’’ cancellation is best understood in a general context [10]. But the upshot for the present purpose is that these higher order eigenfunctions need not be evaluated. 5. Strain energy In terms of the stress components, the strain energy of the beam is given by Ref. [11], article 90: ZZ L 1 n 1 2 ðs 2x þ s 2y þ s 2z Þ ðs x s y þ s y s z þ s z s x Þ þ ðt xy þ t 2yz þ t 2zx Þ dx dy, U¼ 2 2E 2G E
(33)
in which a trivial integration over z has been carried out. Denote U ¼ U0 þ U2 q2 þ :
(34)
5.1. Zeroth order ð0Þ2 This will involve terms such as txy through to szð0Þ2 , as well as terms such as sxð0Þ sð0Þ y , but since all zeroth-order stress components are zero, one immediately has U0 = 0.
5.2. First order ð1Þ ð0Þ ð1Þ ð0Þ ð1Þ ð1Þ ð0Þ This will involve terms such as 2tð0Þ xy txy through to 2sz sz , as well as terms such as sx sy þ sx sy , but since all zeroth-order and first-order stress components are zero, one immediately has U1 =0.
5.3. Second order ð0Þ ð2Þ ð1Þ2 ð0Þ ð2Þ This will involve terms such as tð1Þ2 þ2szð0Þ sð2Þ xy þ2txy txy through to sz z , as well as terms such as sx sy þ ð1Þ ð2Þ ð0Þ s sy þ sx sy ; again all zeroth-order and first-order stress components are zero, and one immediately has U2 =0. ð1Þ x
5.4. Third order ð3Þ ð1Þ ð2Þ ð0Þ ð3Þ ð1Þ ð2Þ ð0Þ ð3Þ This will involve terms such as 2tð0Þ xy txy þ 2txy txy through to 2sz sz þ 2sz sz , as well as terms such as sx sy þ ð2Þ ð2Þ ð1Þ ð3Þ ð0Þ sð1Þ s þ s s þ s s . Again one immediately has U = 0. 3 x y x y x y
3494
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
5.5. Fourth order ð0Þ ð4Þ ð1Þ ð3Þ ð2Þ2 ð4Þ ð1Þ ð3Þ This will involve terms such as tð2Þ2 þ 2sð0Þ xy þ 2txy txy þ 2txy txy through to sz z sz þ 2sz sz , as well as terms such as ð4Þ ð1Þ ð3Þ ð2Þ ð2Þ ð3Þ ð1Þ ð4Þ ð0Þ sð0Þ x sy þ sx sy þ sx sy þ sx sy þ sx sy . The only non-zero contributor is
sð2Þ z ¼ Ex,
(35)
U4 ¼ EIyy L=4:
(36)
and one finds
5.6. Fifth order ð0Þ ð5Þ ð4Þ ð2Þ ð3Þ ð0Þ ð5Þ ð1Þ ð4Þ ð2Þ ð3Þ txy þ2tð1Þ This involves terms such as 2txy xy txy þ 2txy txy through to 2sz sz þ2sz sz þ 2sz sz , as well as terms such ð5Þ ð1Þ ð4Þ ð2Þ ð3Þ ð3Þ ð2Þ ð4Þ ð1Þ ð5Þ ð0Þ as sð0Þ x sy þ sx sy þ sx sy þ sx sy þ sx sy þ sx sy . The only possible contributors come from the second- and thirdð3Þ ð3Þ order stresses, but since there is no product involving sð2Þ z and txz or tyz , one immediately has U5 = 0. This result (and the same for U1 and U3) is anticipated since the strain energy must be even in q.
5.7. Sixth order ð3Þ2 At this level, there are a variety of terms which do contribute; these are tð3Þ2 xz and tyz , which are straightforward, and ð2Þ ð4Þ ð2Þ ð4Þ ð2Þ ð4Þ also 2sð4Þ z sz , sx sz and sy sz . This suggests that one must determine the direct stress components si , i= x, y, z, but in fact one only needs a knowledge of w3. In terms of the displacement components, one has ð4Þ sð4Þ þ2G x ¼ le
qu4 , qx
(37a)
ð4Þ sð4Þ þ2G y ¼ le
qv4 , qy
(37b)
ð4Þ sð4Þ z ¼ le 2Gw3 ,
(37c)
where eð4Þ ¼
qu4 qv4 þ w3 : qx qy
(38)
The relevant expression in the integrand in U is 1 n ð4Þ ð2Þ ð2Þ ð4Þ ð2Þ ð2sð4Þ z sz Þ ðsy sz þ sy sz Þ, 2E E
(39)
sð2Þ z w3 ,
(40)
which reduces to
where w3 comes from the last term in (19c) and all reference to eigenfunctions beyond third order has disappeared. One then finds EL 1 EL J1 þ J, (41) J2 U6 ¼ 2 4ð1 þ nÞ 2 in terms of two integrals with dimensions of (length)6 defined in terms of the Saint-Venant flexure function: ZZ J1 ¼ xðw þ xy2 Þ dx dy, 2 ) ZZ ( qw n 2 n 2 2 qw þ x þ 1 y þð2 þ nÞxy J2 ¼ þ dx dy: qx 2 2 qy
(42a)
(42b)
6. Kinetic energy The Rayleigh quotient Q (see next section) is essentially the ratio of the strain energy U to the kinetic energy T. (This term is used as a shorthand. The kinetic energy is actually o2 T.) Since U starts with q4 and we want Q only to q6, it suffices to calculate T to just q2. Now we have ZZ rL ðu 2 þv 2 þ w 2 Þ dx dy, T¼ (43) 4
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
3495
where again the trivial integration over z has been carried out. The integrand is ðu0 þ u2 q2 þ Þ2 þ ðv0 þv2 q2 þ Þ2 þðw1 qþ w3 q3 Þ2 : The integral of
2 u0
0
provides the q term, while
2 w1 +2u0u2
(44)
2
provides the q term, and one finds
T ¼ T0 þT2 q2 þ ,
(45)
where
T2 ¼
rL 4
rAL
,
(46a)
½Iyy þ nðIxx Iyy Þ:
(46b)
T0 ¼
4
7. Rayleigh quotient It is well known that the eigenvalue o2 is given by the Rayleigh quotient [12]
o2 ¼ Q
U U4 q4 þ U6 q6 þ ¼ , T T0 þT2 q2 þ
(47)
where U and T are evaluated for the corresponding eigenfunction; from this one finds the coefficients defined by (1) as U4 EIyy , ¼ T0 rA
(48a)
U6 U4 T2 E 2 ¼ f2AJIyy ½Iyy þ nðIxx Iyy Þg, T0 rA2 T0
(48b)
A4 ¼
A6 ¼ so that compared with (10b), we find
k¼
2 2ð1þ nÞIyy
2AJ 1 þ AJ 2 =½2ð1 þ nÞ þ nIyy ðIyy Ixx Þ
:
(49)
This expression is similar but not identical to that given by Hutchinson [7]. But using an identity presented in the discussion of Ref. [7], we have J2 ¼ nJ3 2ð1 þ nÞJ1 ,
(50)
where J3 ¼
ZZ 2 2 x y qw n 2 n qw þ x þ 1 y2 þ xy þ ð2 þ nÞxy dx dy: 2 qx 2 2 qy
(51)
Thus we can finally render (49) into
k¼
2 4ð1þ nÞ2 Iyy : 2ð1 þ nÞAJ 1 þ nAJ 3 þ 2nð1 þ nÞIyy ðIyy Ixx Þ
(52)
In this form, k agrees exactly with the expression for the shear coefficient presented in Refs. [4,5], thus proving the latter without having to resort to the physical assumption of TBT. 8. Discussion 8.1. General remarks Our derivation relies on a single approximation, namely that the wavelength is long, i.e., q is small, so that a power series in q makes sense. There is no need to guess, on physical grounds, what degrees of freedom must be kept. In fact, the new variable g in TBT emerges automatically, in the following way. First, the centre line of the beam tilts by an angle (dropping common time-dependent factors from (11)) qu ¼ qucosqz ¼ qðu0 þ u2 q2 þ Þcosqz: qz On the other hand, the cross-sectional plane tilts by an angle qw qw qw1 qw3 3 ¼ cosqz ¼ qþ q þ cosqz: c¼ qx qx qx qx
(53)
(54)
3496
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
The centre line and the cross-sectional plane will, therefore, deviate from orthogonality by an angle qu qw1 qw3 3 ¼ u0 þ q þ u2 þ q þ cosqz: g ¼ cþ qz qx qx
(55)
Our solution, with no assumptions, shows that the first bracket is zero—recovering the key insight and assumption g ¼ 0 in Euler–Bernoulli theory. The next bracket is not zero, so g must be kept for shorter wavelengths, and moreover its effect is captured once w3 is evaluated, without the need to introduce any parameters. In principle, the expansion in q can be continued and is formally exact, although the evaluation of u4, v4, w5, yis impossible in practice—with two important exceptions to be discussed below. But some general features of the expansion can be noted. First, purely on dimensional grounds, each successive term in (1) has an extra power of (qa)2, where a is a typical transverse dimension of the beam, i.e., A2n ðE=ra2 ÞðqaÞ2n . Second, we do not expect TBT to reproduce the next order exactly, i.e., AT8 aA8 , since for shorter wavelengths, the vibration must be described by more than two variables uðz,tÞ, cðz,tÞ at each z — as can be demonstrated in simple cases (see below). These remarks, taken together, imply that TBT, though highly accurate when qa 5 1, cannot be expected to work when qa 1, since the A8q8yterms become important and cannot be reproduced correctly. The discussion in this paper refers to a simply-supported beam, i.e., hinged–hinged end conditions, but this restriction is unnecessary, since other conditions, e.g., guided–guided or guided–hinged, can be regarded as portions of a multi-span hinged–hinged beam [13]. It should also be mentioned that our ‘‘best’’ choice of k is obtained by matching the q6 term in the dispersion relation of the lowest branch (the one without nodes in the cross-sectional plane), and is of course the optimal one for using TBT to describe oscillations of this type—which are the ones most commonly encountered in engineering practice. If one were interested in other types of oscillations, for example the higher branches, other choices may be more appropriate. 8.2. Solvable examples and possible generalizations In two cases, the expansion in q can actually be carried out to very high orders. A brief discussion is given here, principally to illustrate the qualitative remarks above in a precise setting. First consider flexural vibrations in a hypothetical world of two dimensions, say xz. The partial differential equations in Sections 3 and 4 become ordinary differential equations in x. Moreover, on dimensional grounds, the n th-order eigenfunction must go as ðqaÞn ðx=aÞk , with k r n to ensure regularity when a-0; thus it must be a polynomial in x of maximum order n. The differential equation can be cast into algebraic recursion relations for the polynomial coefficients. With these simplifications, the solution can be carried out to many orders. Next consider longitudinal vibrations in a circular cylinder. Using cylindrical coordinates ðr, y,zÞ and factoring out the trivial y-dependence, one again obtains an ordinary differential equation (though of a slightly more complicated form) in r and for the same dimensional reason, the solution is again a polynomial in r, which likewise can be found to very high orders. We have solved both of these cases in powers of q to 20 orders. The coefficients in the polynomials turn out to be rational functions of n, involving powers up to n3n3 ; though complicated, these can be obtained using algebraic software packages. If one seeks only numerical values for a specific value of n, the computation is much simpler. The results, which can also be checked against exact solutions available in these cases [14] confirm the qualitative features discussed in the last subsection. Incidentally, an extension of the two-dimensional problem is a thin plates of thickness h, which may be treated in a manner similar to that given in the present paper, to provide a derivation of the ‘‘best’’ value of k for the Mindlin theory of thin plates [15]. That study will be given elsewhere. 8.3. Using the Rayleigh quotient The present paper makes use of the Rayleigh quotient, which has the advantage [10] of giving the sixth-order eigenvalue A6 from the third-order eigenfunction w3. Two other nice features should be mentioned as well. First, using any approximate eigenfunction to evaluate U and T, and hence Q, guarantees positivity. In contrast, the approximation o2 ¼ A4 q4 þA6 q6 goes negative for qa = O(1). Second, again on physical grounds, one expects o2 =q2 - constant as q-1 (finite phase velocity). This property is also nicely guaranteed for the Rayleigh quotient, since U involves two extra powers of q compared to T. Because of these properties, the Rayleigh quotient is often accurate over a wider range of qa. For longitudinal vibrations in a circular cylinder, this method leads to 5 percent accuracy up to qa 2:5, namely a cylinder with diameter larger than its length—almost a ‘‘disk’’ rather than a ‘‘rod’’. 9. Conclusion In conclusion, the simple and straightforward strategy of expanding in powers of q provides an alternative method to evaluate the shear coefficient k that is systematic and unambiguous. The expansion, formally exact when carried out to all
K.T. Chan et al. / Journal of Sound and Vibration 330 (2011) 3488–3497
3497
orders, also provides a wider perspective to view TBT: (a) with the ‘‘best’’ value for k, it gives the next OððqaÞ2 Þ correction to the classical Euler–Bernoulli theory, but (b) it is itself unlikely to be accurate when qa 1. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
S.P. Timoshenko, On the transverse vibrations of bars of uniform cross-section, Philosophical Magazine 43 (1922) 125–131. A.E.H. Love, A Treatise on the Mathematical Theory of Elasticity, fourth ed., Dover, New York, 1944. N.G. Stephen, On ‘‘A check on the accuracy of Timoshenko’s beam theory’’, Journal of Sound and Vibration 257 (2002) 809–812. N.G. Stephen, Timoshenko’s shear coefficient from a beam subjected to gravity loading, ASME Journal of Applied Mechanics 47 (1980) 121–127. N.G. Stephen, M. Levinson, A second order beam theory, Journal of Sound and Vibration 67 (1979) 293–305. G.R. Cowper, The shear coefficient in Timoshenko’s beam theory, ASME Journal of Applied Mechanics 33 (1966) 335–340. J.R. Hutchinson, Shear coefficients for Timoshenko beam theory, ASME Journal of Applied Mechanics 68 (2001) 87–92 See also Discussion and Closure, ASME Journal of Applied Mechanics 68 (2001) 959–961. J.R. Hutchinson, Shear coefficients for thin-walled Timoshenko beams, Third International Symposium on the Vibrations of Continuous Systems, Jackson Hole, Wyoming, USA, July 2001. J.D. Renton, A check on the accuracy of Timoshenko’s beam theory, Journal of Sound and Vibration 245 (2001) 559–561. K.T. Chan, N.G. Stephen, K. Young, Perturbation theory and the Rayleigh quotient. Journal of Sound and Vibration, 330 (2011) 2073–2078 (Ref: Ms. No. JSVD-10-00949R1). S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, third ed., McGraw-Hill, 1970. S.S. Rao, Mechanical Vibrations, fourth ed., Pearson Prentice Hall, Upper Saddle River, NJ, USA, 2004. N.G. Stephen, The second spectrum of Timoshenko beam theory—further assessment, Journal of Sound and Vibration 292 (2006) 372–389. A.E. Armenakis, D.C. Gazis, G. Herrmann, Free Vibration of Circular Cylindrical Shells, Pergamon, Oxford, 1969. R.D. Mindlin, Influence of rotatory inertia and shear on flexural motions of isotropic elastic plates, ASME Journal of Applied Mechanics 18 (1951) 31–38.