Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 53, 3, pp. 731-743, Warsaw 2015 DOI: 10.15632/jtam-pl.53.3.731 NON-LINEAR VIBRATION OF TIMOSHENKO BEAMS BY FINITE ELEMENT METHOD Jerzy Rakowski, Michał Guminiak Poznań University of Technology, Institute of Structural Engineering, Poznań, Poland e-mail: jerzy.rakowski@put.poznan.pl; michal.guminiak@put.poznan.pl The paper is concerned with free vibrations of geometrically non-linear elastic Timoshen- ko beams with immovable supports. The equations of motion are derived by applying the Hamilton principle. The approximate solutions are based on the negligence of longitudinal inertia forces but inclusion of longitudinal deformations. The Ritz method is used to de- termine non-linear modes and the associated non-linear natural frequencies depending on the vibration amplitude. The beam is discretized into linear elements with independent di- splacement fields. Consideration of the beams divided into the regular mesh enables one to express the equilibrium conditions for an arbitrary large number of elements in form of one difference equation.Owing to this, it is possible to obtain an analytical solution of the dyna- mic problem although it has been formulated by the finite elementmethod. Some numerical results are given to show the effects of vibration amplitude, shear deformation, thickness ratio, rotary inertia, mass distribution and boundary conditions on the non-linear natural frequencies of discrete Timoshenko beams. Keywords: Timoshenko beams, non-linear vibrations, eigenvalue problem, finite elements 1. Introduction The linear free vibration of Timoshenko beamswas studied byLevinson andCooke (1982). Bha- shyam and Prathap (1981) applied the Finite Element Method and elements with linear shape functions to determine natural frequencies of Timoshenko beams. The authors eliminated the shear locking by a selective reduced integrations procedure. The studies of beamswith immova- ble supports weere carried out by Woinowsky-Krieger (1950), Hsu (1960), Evensen (1968) and Lewandowski (1987). The authors used a continuum approach as well as the finite element ap- proach proposed byLevinson andCooke (1982), Sarma andVaradan (1983) andKitipornchai et al. (2009). The exact solutions given in form of elliptic functions were presented byWoinowsky- -Krieger (1950) and Hsu (1960). Evensen (1968) analysed beams with various boundary condi- tions using the perturbation method. Lewandowski (1987) applied the Ritz method to obtain the frequency-amplitude relationship for elastically supportedbeams.The solutions presentedby Bhashyam and Prathap (1980) as well as Sarma and Varadan (1983) are based on the assump- tion that the equation of motion is satisfied only at the instant of the maximum deformation. Dumir and Bhaskar (1988) compared the results of analysis of non-linear vibrations of beams andplates obtainedbyvariousmethods.Theyalso brought out the source of errorswhich appear in some finite elements formulations, e.g. presented by Mei and Decha-Umphai (1985). Marur and Prathap (2005) solved the non-linear vibration problem of Timoshenko beams by applying variational formulations, and similar isssues for thick asymetric beams were inastigated in a wide range by Singh and Rao (1998). The problem of non-linear vibration of Timoshenko be- ams taking into account cracking was invastigated byKitipornchai et al. (2009). The non-linear formulation of the Timoshenko beam based on themodified couple stress theorywas applied by Asghari et al. (2010) to static and free vibration analysis of beams. 732 J. Rakowski,M. Guminiak Effectiveness of the Finite ElementMethod depends e.g. on the type of finite elements taken in the analysis. Many of them induce undesirable parasitic effects as shear locking. A critical analysis of a class of finite elements was carried out by Rakowski (1990, 1991a, 1991b). In the present work the linear finite element developed by Rakowski (1990) is used in the analysis of non-linear vibration of theTimoshenkobeam.Suchan improved linear element doesnot lock and gives a better convergence to the exact results than the reduced-integrated one, and the division of considered beams into the regular mesh of identical elements enables one to formulate the equilibriumconditions in formofonedifference equation.Theadvantage of thisFEMformulation idea is that the solution of the considered dynamic problem can be obtained in an analytical closed form and that it facilitates analyses regardless of the number of finite elements. The solving procedure is analogous to that used for continuous systems, but in this case the nodal displacements are spatial functions of the discrete variable (the index number of a node). This methodology was adopted to solve the problem of vibrations of infinite Bernoulli-Euler beams by Rakowski andWielentejczyk (1996, 2002). In the present approach, the equations of motion for a continuous beam are derived from the variational formulation by applying Hamilton’s principle. The approximate time function is assumed to be harmonic and the solution of the problem is based on the idea presented by Rosenberg (1966) and Szemplińska-Stupnicka (1983) wherenon-linearnormalmodesandnatural frequencies of vibrations areunknownandamplitude dependent. The effects of shear deformation, thickness ratio and rotary inertia on the non-linear natural frequencies are analysed. The element stiffness matrix and the consistent mass matrix are derived using the simple linear elements with independent displacement fields. The shear locking in exactly integrated elements is eliminated by introducing a scalling factor into the stiffness matrix (Rakowski, 1990). Beams with different boundary conditions are considered. The comparison of results obtained for the evaluated elements and the reduced-integrated ones is given. 2. Theoretical considerations The equation for the strain energy for a beam of length l, including shear deformation, axial force, longitudal and the large transverse displacement effect, is U = l ∫ 0 EA 2 ( u,x− 1 2 w2,x )2 dx+ l ∫ 0 EJ 2 θ2,x dx+ l ∫ 0 κGA 2 (w,x−θ) 2 dx + l ∫ 0 H ( u,x− 1 2 w2,x )2 dx (2.1) where u, w are the longitudinal and transverse displacements of the centroidal axis, θ is the rotation of the beam section, E, G are the elastic and shear moduli, J, A are the moment of inertia and the cross-section area, respectively, H is the initial axis force and κ is the shear coefficient. The kinetic energy including horizontal and rotary inertias can be expressed as T = l ∫ 0 m 2 u2,t dx+ l ∫ 0 m 2 w2,t dx+ l ∫ 0 mρ2 2 θ2,t dx (2.2) wherem is mass per unit length and ρ is the radius of gyration of the cross section. Equations (2.1) and (2.2) give, on applying Hamilton’s principle, the following equations of motion Non-linear vibration of Timoshenko beams by FEM 733 EA [ u,x+ 1 2 (w,x) 2 ] ,x −mu,tt =0 { Hw,x+EA [ u,x+ 1 2 (w,x) 2 ] w,x } ,x +κGA[w,x−θ],x−m,tt =0 EJθ,xx+κGA(w,x−θ)−mρ 2θ,tt =0 (2.3) Neglecting the horizontal inertia forces in the second equation of (2.3), yields C =EA [ u,x+ 1 2 (w,x) 2 ] ,x = const (2.4) Integrating expression (2.4) in the range from 0 to l for an immovable beam, we have C = l ∫ 0 EA 2a w2,x dx (2.5) Substituting (2.4) into equation (2.3) and eliminating the unknown function θ, yields (H+C κGA +1 )∂4w ∂x4 − H+C EJ ∂2w ∂x2 − [ m κGA + mρ2 EJ (H+C κGA +1 )] ∂4w ∂x2∂t2 + m EJ ∂2w ∂t2 + m2ρ2 κEJGA ∂4w ∂t4 =0 (2.6) Let us assume the approximate solution to equation (2.6) to be harmonic function of time (Woinowsky-Krieger, 1950) w(x,t) =αW(x)cosωt (2.7) whereW(x) is a function of the spatial variable (not given a priori),α is the vibration amplitude of anarbitrarypointx0 of thebeamaxis forwhichW(x0)= 1.Thenon-linearnormalmodeW(x) and the free frequency vibrations ω depend on the boundary conditions and on the vibration amplitude α. The above solution does not satisfy non-linear equation (2.6) at every time. We minimize the error by equating theweighed residual of equation (2.6) to zero for one time period with theweighting function cosωt (Dumir andBhaskar, 1988; Szemplińska-Stupnicka, 1983). As a result, the equilibrium equation in the dimensionless form is obtained d4W dξ4 + 1 F [ Ω2 F +f λ2 − (s+ c)λ2 ]d2W dξ2 + Ω2 F ( f Ω2 λ4 −1 ) W =0 (2.8) where W =W(ξ) ξ= x L f = 2(1+ν) κ s= H EA λ= L ρ Ω= mω2L4 EJ c= 3 8 δ2 λ2 1 ∫ 0 dW dξ dξ δ= α ρ F = s+ c f +1 and ν is Poisson’s ratio. The solution tf equation (2.8) is W(x)=C1 sinhµ1ξ+C2coshµ1ξ+C3 sinµ2ξ+C4cosµ2ξ (2.9) whereCi are arbitrary constants and µ21 = √ b2 4 − Ω F (fΩ2 λ4 −1 ) + b 2 µ22 = √ b2 4 − Ω F (fΩ2 λ4 −1 ) − b 2 b= 1 F [Ω2 λ2 (F +f)− (s+ c)λ2 ] 734 J. Rakowski,M. Guminiak From the third equation of (2.3) one finds that thedimensionless rotation functiondepending on the spatial variable is of the form φ(ξ)= (C1 sinhµ1ξ+C2coshµ1ξ) β1 µ1 − (C3cosµ2ξ+C4 sinµ2ξ) β2 µ2 (2.10) where β1 = fΩ 2/λ2+Fµ21 and β2 = fΩ 2/λ2−Fµ22. The boundary conditions for beams with immovable supports are: — for the hinged end W =0 M =EJ dφ dξ =0 — for the clamped end W =0 φ=0 The frequency equations derived for several combinations of the boundary conditions are: — for the hinged-hinged beam sinhµ1 sinµ2 =0 (2.11) — for the clamped-clamped beam 2+ β22µ 2 1−β 2 1µ 2 2 β1β2µ1µ2 sinhµ1 sinµ2−2coshµ1cosµ2 =0 (2.12) — for the clamped-hinged beam β2µ1 tanhµ1+β1µ2 tanhµ2 =0 (2.13) Since the non-linear frequency parameter Ω2 and the non-linear vibration modeW depend on the vibration amplitude α, the eigenvalue problem of non-linear equation (2.8) is solved iteratively (Lewandowski, 1987). For hinged-hinged beams, the non-linear mode is identical to the linear one and independent of the amplitude. Hence the frequency parameter Ω2 can be expressed in an analytical form as the root of the quadratic equation f Ω4 λ4 −Ω2 [ (F +k) k2π2 λ2 +1 ] +k2π2[Fk2π2+λ(s+ c)] = 0 (2.14) where c = (3/16)(k2π2δ2/λ2), k is the number of the vibration mode. Ommiting the rotary inertia effect, yields Ω2 = k2π2 Fk2π2+λ2(s+ c) 1+f k 2π2 λ2 (2.15) By assuming in (2.14) or (2.15) f =0, one obtains non-linear frequencies for the Bernoulli- -Euler beam in which the axial force effect is taken into account. Table 1 presents the rotary inertia effect Iρ on the non-linear frequency ratio ω/ω1 for the fundamental vibration mode of simply-supported beams (α is the amplitude in the middle of the beam, Ω1 is the frequency parameter determined for linear vibration). The calculations have been carried out for ν =0.3, κ=5/6 andH =0, respectively. The backbone curves for various slenderness ratios for two modes of vibration of a simply- -supported Timoshenko beam are shown in Fig. 1. For the firstmode (k=1),α is the amplitude in themiddle of the beam (ξ=1/2), for k=2, ξ =1/4 (L is the beam length). One can see that the frequency ratio for the simply-supported Timoshenko beam depends on the kind of mode. Non-linear vibration of Timoshenko beams by FEM 735 Table 1.Non-linear frequency ratio ω/ω1 for simply-supported continuous beams Bernoulli-Euler beam Timoshenko beam α/ρ λ=30 λ=10 λ=20 λ=30 Iρ 6=0 Iρ =0 Iρ 6=0 Iρ =0 Iρ 6=0 Iρ =0 Iρ 6=0 Iρ =0 1.0 1.0897 1.0897 1.1158 1.1159 1.0963 1.0963 1.0927 1.0927 2.0 1.3229 1.3229 1.4068 1.4075 1.3445 1.3445 1.3325 1.3325 3.0 1.6394 1.6394 1.7889 1.7908 1.6785 1.6785 1.6559 1.6559 4.0 2.0000 2.0000 2.2146 2.2190 2.0568 2.0569 2.0255 2.0255 Ω1 9.8159 9.8696 8.3874 8.6299 9.4106 9.5103 9.6556 9.7050 Fig. 1. The backbone curves for various slenderness ratios for twomodes of vibration of a simply-supported Timoshenko beam 3. Finite Element Method formulation Let us assume the element shape functions to be linear and independent u(x,t) =u1(t)+ [u2(t)−u1(t)] x a w(x,t) =w1(t)+ [w2(t)−w1(t)] x a θ(x,t)= θ1(t)+ [θ2(t)−θ1(t)] x a (3.1) The beam element convention is presented in Fig. 2. Substituting (3.1) into (2.1), we obtain the element stiffness matrix in which the non-linear part corresponding to the nodal quantities {u1 u2 w1 θ1 w2 θ2} T can be written in the following form 736 J. Rakowski,M. Guminiak Fig. 2. Timoshenko beam convention K ∗ N= EA a          1 −1 (w2−w1)/(2a) 0 −(w2−w1)/(2a) 0 −1 1 −(w2−w1)/(2a) 0 (w2−w1)/a 0 (w2−w1)/a −(w2−w1)/a (w2−w1) 2/(2a2) 0 −(w2−w1) 2/(2a2) 0 0 0 0 0 0 0 −(w2−w1)/a (w2−w1)/a −(w2−w1) 2/(2a2) 0 (w2−w1) 2/(2a2) 0 0 0 0 0 0 0          Neglecting the horizontal forces in the calculations yields for each element the following rela- tionship C = EA a [ u1−u2+ w1(w2−w1) 2a − w2(w2−w1) 2a ] = const (3.2) Assuming that the considered Timoshenko beam with immovable supports is divided into elements with equally spaced nodes and summing up both sides of equation (3.2) in the range from 0 toR (R is the number of beam elements), we have C = EA 2Ra R−1 ∑ r=0 (wr+1−wr) 2 (3.3) In order to avoid the shear lockingwhich occurs in linear elements and leads to erroneous results, we improve the elements by introducing a scaling factor d/(d+1) into the flexural and shear matrices (Rakowski, 1990). ReplacingE andG byEd/(d+1) andGd/(d+1), respectively, and eliminating the unknowns ui, the element stiffness matrix can be expressed as follows K= EJ a3 1 d+1 [dKF +KS +(d+1)(KG+KN)] (3.4) where d = [24(1+ν)ρ2]/(κa4), KF , KS, KG and KN are the flexural, shear, geometrical and non-linear matrices, respectively. They refer to the nodal quantities {w1 aθ1 w2 aθ2} T and have the form KF =      0 0 0 0 0 1 0 −1 0 0 0 0 0 −1 0 1      KS =      12 6 −12 6 6 4 −6 2 −12 −6 12 −6 6 2 −6 4      KG = a2 EJ      H 0 −H 0 0 0 0 0 −H 0 H 0 0 0 0 0      KN = a2 EJ      C 0 −C 0 0 0 0 0 −C 0 C 0 0 0 0 0      (3.5) Non-linear vibration of Timoshenko beams by FEM 737 We assume the nodal displacements to be harmonic wr =αWr cosωt θr = α a φr cosωt (3.6) whereα is the amplitude of an arbitrarly chosen node i,Wr and φr are dimensionless functions of the discrete variable r (r is the index number of the node), r ∈ 〈0,R〉. By using the shape functions given in (3.1), the consistent matrix is obtained as MC = am 6      2 0 1 0 0 2ρ2/a2 0 ρ2/a2 1 0 2 0 0 ρ2/a2 0 2ρ2/a2      (3.7) If the beam, divided into a regular mesh, consists of identical elements, the dynamic equili- brium conditions can be derived according to the following way. We substitute (3.5) into (3.4) and take into account element mass matrix (3.7). Having assembled the two adjacent beam elements (r− 1,r) and (r,r+1), we can write the equilibrium of moments and forces for an arbitrary node r as 6(Wr−1−Wr)+2(φr−1+2φr)+6(Wr −Wr+1)+2(2φr +φr+1)−d(φr−1+φr) +d(φr −φr+1)−4Λ 2(d+1)e2(φr−1+4φr +φr+1)=0 − [12+6(d+1)](Wr−1−Wr)+ [12+6(d+1)D](Wr −Wr+1)−6(φr−1−φr) +6(φr −φr+1)−4Λ 2(d+1)(Wr−1+4Wr+Wr+1)= 0 (3.8) Introducing the shifting operator Enr and central difference operator∆ n r into equations (3.8), we obtain two difference equations with unknown nodal variables φr andWr [ (∆2+6)− d 2 ∆2 ] φr−3(E−E −1)Wr −2Λ 2(d+1)e2(∆2+6)φr =0 (E−E−1)φr−2[2+D(d+1)]∆ 2Wr− 2 3 Λ2(d+1)(∆2+6)Wr =0 (3.9) where ∆2 =∆2r =E+E −1−2 En =Enr Λ= Ω2 24R4 Ω2 = mω2L4 EJ D= a2(H+C) 6EJ e= R λ L is length of the beam (L= aR), λ is the slenderness ratio (λ=L/ρ. Equating to zero the weighed residual of equation (3.9) for one time period (0− 2π) with the weighting function cosωt (Marur and Prathap, 2005) 2π ∫ 0 εi(r,t)cosωt d(ωt)= 0 i=1,2 (3.10) (ε1, ε2 are the left-hand sides of the first and second equation of (3.9) multiplied by cosωt) one obtains the equilibrium conditions given in the form of equations (3.9) where now D= a2 6EJ ( H+ 3 4 C ) C = EARδ2 2λ2 R−1 ∑ r=0 (Wr+1−Wr) 2 (3.11) 738 J. Rakowski,M. Guminiak After elimination of φr from coupled equation (3.9), we can rewrite them in form of a single fourth-order homogenous difference equation B1∆ 4Wr+B2∆ 2Wr+B3Wr =0 (3.12) where B1 =1−D ( 1− d 2 ) + Λ2 3 (d−2+12e2)+2Λ2De2(d+1)+ 4 3 Λ4e2(d+1) B2 =−6D+2Λ 2(d−4+12e2)+12Λ2De2(d+1)+16Λ4e(d+1) B3 =−24Λ 2+48Λ4e2(d+1) The analytical solution to equation (3.12) is Wr =C1 sinhµ1r+C2coshµ1r+C3 sinhµ2r+C4coshµ2r (3.13) whereµ1,µ2 can be real, imaginary or complex, and theymust fulfill the characteristic quadratic equation B1cosh 2µ+ ( −2B1+ B2 2 ) coshµ+B1−B2+ B4 4 =0 (3.14) The function φr can be derived from the second equation of (3.9) in the form φr =D2 sinhµ1r+D1coshµ1r+D4 sinhµ2r+D3coshµ2r (3.15) where D1 =β1C1 D2 =β1C2 D3 =β2C3 D4 =β2C4 β1 = 2 sinhµ1 { (coshµ1−1) [ 1+ D 2 (d+1) ] + Λ2 3 (coshµ1+2)(d+1) } β2 = 2 sinhµ2 { (coshµ2−1) [ 1+ D 2 (d+1) ] + Λ2 3 (coshµ2+2)(d+1) } If a lumpedmassmodel of theTimoshenko beam is considered, the elementmassmatrix has the following form M1 = am 2      1 0 1 0 0 ρ2/+a2 0 ρ2/a2 1 0 1 0 0 ρ2/a2 0 ρ2/a2a2      (3.16) The difference equilibrium equations corresponding to element mass matrix M1 (3.16) are [ (∆2+6)− d 2 ∆2 ] φr−3(E−E −1)Wr−12Λ 2e2(d+1)φr =0 (E−E−1)φr −2[2+D(d+1)]∆ 2Wr−4∆ 2(d+1)Wr =0 (3.17) or expressed in form of an equation with one unknownWr B1∆ 4Wr+B2∆ 2Wr+B3Wr =0 where B1 =1−D ( 1− d 2 ) B2 =−6D+2Λ 2[d−2+12e2+6De2(d+1)] B3 =−24Λ 2+48Λ4e2(d+1) Non-linear vibration of Timoshenko beams by FEM 739 The solution to this equation and the expression of the characteristic equation are identical to those given in (3.13) and (3.14), respectively. But in this case, the quantities βi present in the relationship between Di andCi (see (3.15)), derived from (3.17), are defined as β1 = 2 sinhµi { coshµi−1) [ 1+ D 2 (d+1) ] +Λ2(d+1) } i=1,2 Since the non-linear frequency parameter Λ and the non-linear vibration modes Wr and φr depend on the amplitude α, the eigenvalue problem of non-linear equation (3.12) is solved ite- ratively. The iteration starts by assuming the magnitude of the constant D (it is convenient in the first step of iteration to use the value calculated for the linear case). From the frequency equations adequate to the considered boundary conditions, one finds µ1 and µ2. Then the vi- bration mode Wr is determined and the resulting parameter D is compared with the given D. The iterative computation is continued until the difference |D−D| is less than the adopted tolerance. The boundary conditions are considered below. 3.1. Simply-supported beam For a simply supported beam with the boundary nodes 0 and R, the boundary conditions W0 =WR =0,M0 =MR =0 lead to the following mode of vibrations which is independent of the vibration amplitude Wr =C1 sin(kπr) (3.18) and k being the number of mode. The frequency parameterΛ can be determined as the root of the quadratic equation akΛ 2+ bkΛ+ ck =0 (3.19) where ak =4e 2(d+1) {1 3 [ cos (kπ R ) −1 ]2 +2 [ cos (kπ R ) −1 ] +3 } bk = [ cos (kπ R ) −1 ]2[ − 2 3 + d 3 +4e2+2e2(d+1)D ] + [ cos (kπ R ) −1 ] [−4+d+12e2+6e2(d+1)D]−6 ck = [ cos (kπ R ) −1 ]{[ cos (kπ R ) −1 ][ 1−D ( 1− d 2 )] −3D } When the lumpedmass model is considered, the above coefficients must be replaced by ak =4e 2(d+1) bk = [ cos (kπ R ) −1 ]2 [−2+d+12e2+6e2(d+1)D]−6 ck = [ cos (kπ R ) −1 ]{[ cos (kπ R ) −1 ][ 1−D ( 1− d 2 )] −3D } 3.2. Clamped-clamped beam The boundary conditions W0 =WR =0, φ0 =φR =0 request the following equations to be fulfilled C1[β2 sinh(µ1R)−β1 sinh(µ2R)]+C2β2[cosh(µ1R)− cosh(µ2R)] = 0 C1β1[cosh(µ1R)− cosh(µ2R)]+C2[β1 sinh(µ1R)−β2 sinh(µ2R)] = 0 (3.20) 740 J. Rakowski,M. Guminiak Equating to zero the determinant of the set of equations (3.20) yields the characteristic equation 2β1β2[cosh(µ1R)cosh(µ2R)−1]− (β 2 1 +β 2 2)sinh(µ1R)sinh(µ1R)= 0 (3.21) Themodes of vibrations are Wr =C1 { sinh(µ1r)− β1 β2 sinh(µ2r)−β[cosh(µ1r)− cosh(µ2r)] } (3.22) where β= β1[cosh(µ1R)− cosh(µ2R)] β1 sinh(µ1R)−β2 sinh(µ2R) 3.3. Hinged-clamped beam In this case, the boundary conditions are given by W0 = WR = 0, M0 = 0, φR = 0, which refer to the homogenous equations C1 sinh(µ1R)+C3 sinh(µ2R)=0 C1β1cosh(µ1R)+C3β2cosh(µ2R)= 0 (3.23) and give the characteristic equation β2 tanh(µ1R)−β1 tanh(µ2R)= 0 (3.24) Themodes of vibrations can be expressed as Wr =C1 [ sinh(µ1r)− sinh(µ1R) sinh(µ2R) sinh(µ2r) ] (3.25) Let us now analyse the reduced-integrated linear elements. When the one-point quadrature is applied to the shear term of strain energy (2.1), the element stiffness matrix becomes K= EJ a3 1 d [(d−1)KF +KS +d(KG+KN)] (3.26) where KF , KS, KG and KN are defined as in (3.4). The difference between the above matrix and that given in (3.4) is that in (3.26) the expression (d−1) is used instead of d. 4. Numerical examples Thenon-linear frequency ratios are determined for beamswith various boundary conditions and for a wide spectrum of the thickness ratio. The consistent and lumped mass model of beams are considered. The results of computation are compared with the results obtained for linear elements towhich the reduced-selective integration technique has been applied.The calculations have been carried out with the values ν = 0.3, κ = 5/6 and H = 0. The following indications are accepted: • CONT – the beamwith continuous mass distribution, • MOD – the beam divided into modified finite elements proposed by Rakowski (1990), • RI – the beam divided into linear elements with reduced integration correction. Numerical results are obtained for the lumped mass model. Values in the parentheses refer to the consistent mass model. Non-linear vibration of Timoshenko beams by FEM 741 Table 2.Comparison of the frequency ratio (ω/ω1) for continuous and discrete (R=8) simply- -supported Timoshenko beams α/ρ λ=20 λ=100 CONT MOD RI CONT MOD RI First mode (k=1) Ω1 9.411 9.528 9.584 9.850 9.977 10.042 (9.406) (9.461) (9.849) (9.914) 1.0 1.0963 1.0941 1.0930 1.0900 1.0878 1.0867 2.0 1.3445 1.3371 1.3336 1.3237 1.3165 1.3129 3.0 1.6785 1.6652 1.6587 1.6409 1.6278 1.6212 4.0 2.0568 2.0375 2.0281 2.0023 1.9832 1.9736 Secondmode (k=2) Ω1 33.550 35.115 35.772 39.162 41.201 42.338 (33.357) (33.981) (39.138) (40.218) 1.0 1.1158 1.1062 1.1021 1.0908 1.0823 1.0781 2.0 1.4048 1.3763 1.3629 1.3264 1.2983 1.2841 3.0 1.7889 1.7346 1.7108 1.6457 1.5945 1.5685 4.0 2.2146 2.1363 2.1021 2.0092 1.9345 1.8962 Table 3. Frequency ratio (ω/ω1) of the first mode for the clamped-clamped beam (R=8) α/ρ λ=20 λ=40 λ=100 MOD RI MOD RI MOD RI 1.0 1.0251 1.0242 1.0211 1.0204 1.0202 1.0195 (1.0251) (1.0242) (1.0212) (1.0294) (1.0202) (1.0195) 2.0 1.0963 1.0929 1.0817 1.0789 1.0781 1.0756 (1.0965) (1.0931) (1.0819) (1.0790) (1.0783) (1.0757) 3.0 1.2045 1.1975 1.1746 1.1689 1.1675 1.1625 (1.2049) (1.1979) (1.1750) (1.1693) (1.1678) (1.1628) 4.0 1.3400 1.3286 1.2918 1.2830 1.2807 1.2734 (1.3407) (1.3294) (1.2925) (1.2836) (1.2814) (1.2739) 5.0 1.4952 1.4789 1.4265 1.4144 1.4114 1.4021 (1.4962) (1.4800) (1.4276) (1.4154) (1.4123) (1.4028) Ω1 19.082 19.539 21.609 22.281 22.529 23.295 (18.801) (19.247) (21.282) (21.937) (22.185) (22.932) CONT Ω1 18.837 21.296 22.189 5. Concluding remarks The analytical solution to the non-linear free vibration problem formulated by the finite element method for theTimoshenko beam is presented.The simplemodified linear element is used in the solution. This improved element gives a better convergence to the solution for the continuous system than the reduced-integrated one. The good accuracy is achieved for a wide range of the thickness ratio for the consistent and lumped mass models of the beam. The concept of the finite difference method is introduced. The main advantage of the presented idea is that even for a large number of elements it suffices to consider the non-linear eigenvalue problem of only one difference equation which is equivalent to a set of FEM equilibrium conditions. 742 J. Rakowski,M. Guminiak Table 4. Frequency ratio (ω/ω1) of the secondmode for the clamped-clamped beam (R=8) α/ρ λ=20 λ=40 λ=100 MOD RI MOD RI MOD RI 1.0 1.0649 1.0604 1.0452 1.0413 1.0402 1.0367 (1.0652) (1.0608) (1.0457) (1.0418) (1.407) (1.10371) 2.0 1.2372 1.2216 1.1680 1.1540 1.1502 1.1379 (1.2383) (1.2228) (1.1697) (1.1557) (1.1519) (1.1393) 3.0 1.4780 1.4480 1.3431 1.3154 1.3079 1.2845 (1.4796) (1.4499) (1.3462) (1.3188) (1.3114) (1.2875) 4.0 1.7591 1.7134 1.5506 1.5068 1.4950 1.4590 (1.7611) (1.7158) (1.5552) (1.5120) (1.5004) (1.4641) 5.0 2.0640 2.0024 1.7790 1.7169 1.7003 1.6499 (2.0662) (2.0050) (1.7846) (1.7236) (1.7074) (1.6572) Ω1 46.229 48.146 58.332 62.197 63.945 69.172 (43.863) (45.585) (55.152) (58.717) (60.409) (65.225) CONT Ω1 44.330 55.423 60.519 Table 5. Frequency ratio (ω/ω1) of the third mode for the clamped-clamped beam (R=8) α/ρ λ=20 λ=40 λ=100 MOD RI MOD RI MOD RI 1.0 1.1095 1.0989 1.0586 1.0473 1.0442 1.0335 (1.1117) (1.1009) (1.0603) (1.0489) (1.0457) (1.0347) 2.0 1.3768 1.3433 1.2124 1.1737 1.1630 1.1254 (1.3799) (1.3461) (1.2126) (1.1777) (1.1671) (1.1290) 3.0 1.7252 1.6652 1.4268 1.3532 1.3329 1.2599 (1.7262) (1.6657) (1.4310) (1.3581) (1.3384) (1.2652) 4.0 2.1172 2.0290 1.6793 1.5673 1.5367 1.4234 (2.1133) (2.0241) (1.6825) (1.5717) (1.5425) (1.4300) 5.0 2.5518 2.4235 1.9562 1.8047 1.7774 1.6234 (2.5283) (2.4042) (1.9579) (1.8082) (1.7694) (1.6152) Ω1 81.697 85.847 112.487 124.642 130.900 151.802 (72.825) (76.475) (99.927) (110.404) (116.090) (133.961) CONT Ω1 75.077 101.659 117.002 References 1. Asghari M., Kahrobaiyan M.H., Ahmadian M.T., 2010, A nonlinear Timoshenko beam for- mulation based on themodified couple stress theory, International Journal of Engineering Science, 48, 1749-1776 2. Bhashyam G.R., Prathap G., 1980, Galerkin finite elementmethod for non-linear beam vibra- tions, Journal of Sound and Vibration, 72, 191-203 3. Bhashyam G.R., Prathap G., 1981, The second frequency spectrum of Timoshenko beams, Journal of Sound and Vibration, 76, 407-420 4. Dumir P.C., Bhaskar A., 1988, Some erroneous finite element formulations of non-linear vibra- tions of beams and plates, Journal of Sound and Vibration, 123, 517-527 Non-linear vibration of Timoshenko beams by FEM 743 5. EvensenD.A., 1968,Non-linear vibrations of beamswith various boundary conditions,American Institute of Aeronatics and Astronomics Journal, 6, 370-372 6. HsuC.S., 1960,On the application of elliptic functions in non-linear forced oscillations,Quarterly of Applied Mathematics, 17, 393-407 7. Kitipornchai S., Ke L.L., Yang J., Xiang Y., 2009, Nonlinear vibration of edge cracked functionally graded Timoshenko beams, Journal of Sound and Vibration, 324, 962-982 8. LewandowskiR., 1987,Applicationof theRitzmethod to the analysis of non-linear free vibration of beams, Journal of Sound and Vibration, 114, 91-101 9. Levinson M., Cooke D.W., 1982, On the two frequency spectra of Timoshenko beams, Journal of Sound and Vibration, 84, 319-326 10. MarurS.R., PrathapG., 2005,Non-linearbeamvibrationproblems and simplifications in finite element models,Computational Mechanics, 35, 352-360 11. Mei C., Decha-Umphai K., 1985, A finite element method for non-linear forced vibrations of beams, Journal of Sound and Vibration, 102, 369-380 12. Rakowski J., 1990, The interpretation of shear locking in beam elements,Computers and Struc- tures, 37, 5, 769-776 13. Rakowski J., 1991a, A critical analysis of quadratic beam finite elements, International Journal for Numerical Methods in Engineering, 31, 949-966 14. Rakowski J., 1991b, A new methodology of evaluation of C0 bending finite elements,Computer Methods in Applied Mechanics and Engineering, 91, 1327-1338 15. Rakowski J., Wielentejczyk P., 1996, Vibrations of infinite periodic beams by finite element method, Zeitschrift für Angewandte Mathemetic und Mechanik, 76, 411-412 16. Rakowski J., Wielentejczyk P., 2002, Dynamic analysis of infinite discrete structures,Foun- dations of Civil and Environmental Engineering, 3, 91-106 17. Rosenberg R.M., 1966, On non-linear vibrations of systems with many degree of freedom, Ad- vanced in Applied Mechanics, 9, 154-242 18. Sarma B., Varadan T.K., 1983, Lagrange-type formulation for finite element analysis of non- linear beam vibrations, Journal of Sound and Vibration, 86, 61-70 19. Singh G., Rao G.V., 1998, Nonlinear oscillations of thick asymmetric cross-ply beams, Acta Mechanics, 127, 135-146 20. Szemplińska-StupnickaW., 1983,Non-linear normalmodes andgeneralizedRitzmethod in the problems of vibrations of non-linear elastic continuous systems, International Journal of Non-linear Mechanics, 18, 154-242 21. Woinowsky-Krieger S., 1950, The effect of an axial force on the vibration of hinged bars, Journal of Applied Mechanics, 17, 35-36 Manuscript received January 30, 2014; accepted for print March 12, 2015