An Investigation on the Nonlinear Free Vibration Analysis of Beams with Simply Supported Boundary Conditions Using Four Engineering Theories

The objective of this study is to present a brief survey on the geometrically nonlinear free vibrations of the Bernoulli-Euler, the Rayleigh, shear, and the Timoshenko beams with simple end conditions using the Homotopy Analysis Method (cid:2) HAM (cid:3) . Expressions for the natural frequencies, the transverse deﬂection, postbuckling load-deﬂection relation to, and critical buckling load are presented. The results of nonlinear analysis are validated with the published results, and excellent agreement is observed. The e ﬀ ects of some parameters, such as slender ratio, the rotary inertia, and the shear deformation, are examined as other parameters are ﬁxed.


Introduction
Many structures, such as bridges, buildings, and spacecraft arms can be modeled as flexible beams. Most phenomena in these structures are essentially nonlinear and are described by nonlinear equations. Therefore, the study of nonlinear effects on vibration analysis of beam structures is of great importance for engineers working in the field of aerospace, mechanical, and civil engineering.
Although the nonlinear free vibrations of classical Bernoulli-Euler's beams have been investigated by many researchers 1-5 , only a few publications have been devoted to include the effects of shear deformation and rotary inertia 6, 7 .
Azrar et al. 1 have developed a semianalytical approach to the nonlinear dynamic response problem based on Lagrange's principle and the harmonic balance method. An analytical method for determining the vibration modes of geometrically nonlinear beams under various edge conditions has been presented by Qaisi 2 . Nonlinear modal analysis approach based on invariant manifold method is utilized to obtain the nonlinear normal modes of a clamped-clamped beam for large amplitude displacements by Xie et al. 3 . Singh et al. 4 have developed a spline-based differential quadrature method based on sextic cardinal spline functions. Guo and Zhong 5 have investigated nonlinear vibrations of thin beams based on sextic cardinal spline functions, a spline-based differential quadrature method. Rao et al. 6 have studied the large-amplitude free vibrations of slender beams using the continuum and finite element methods. The large-amplitude free-vibration behavior of the beams with central point-concentrated masses has investigated by Rao et al. 7 .
The main purpose of this survey is to present the analytical expressions for geometrically nonlinear vibration of Bernoulli-Euler's, Rayleigh's, shear, and Timoshenko's beams with simple end conditions. First it is obtained nonlinear differential equations of motion on the basis of the four engineering beam theories. It is assumed that only the fundamental mode is excited and the governing nonlinear partial differential equation is reduced to a single nonlinear ordinary differential equation. These equations have been solved analytically in time domain using HAM. Finally, some numerical examples are shown for the nonslender and slender beams to signify the differences among the beam models.

Timoshenko's Beam Model
Consider a straight Timoshenko's beam of length L, a uniform cross-sectional area A, the mass per unit length of m, the moment of area of the cross-section I, Young's modulus E, and shear modulus G that is subjected to an axial force of magnitude N as shown in Figure 1. We assume that the beam is made of a homogenous and isotropic material.
The kinetic energy T and the potential energy U of the vibrating beam can be written as Journal of Applied Mathematics 3 where γ represents the shear angle γ w , x − ψ , u, w, and ψ are axial, transverse displacements, and the cross-section rotation due to the bending moment, x is the axial coordinate of the beam, r is the radius of gyration I/A , and k is beam's cross-sectional shape factor. Also M, V , and N x are the bending moment, shear, and axial forces, respectively. Comma denotes differentiate with respect to x or t.
The strain-displacement relation of the Von-Karman type is used: Applying Hamilton's principle, the governing equations of motion and boundary conditions are obtained as follows: Upon neglecting the axial inertia, 2.3 -2.4 can be combined to get the first type of differential equations which is coupled as follows: where the axial force N x is given by In the next step, 2.7 -2.8 can be combined to arrive at the following decoupled equation of motion:

2.10
It may be noted that the equation of motion governing the rotation of the cross-section ψ can be obtained using the same procedures to derive an equation similar to 2.10 with w being replaced by ψ.

Shear Beam Model
This model neglects the effect of rotary inertia in Timoshenko's beam model. Therefore, the kinetic energy can only be rewritten as Applying the same process mentioned in Timoshenko's beam model, the decoupled equation of motion is given by and the boundary conditions are as listed in 2.6 .

Rayleigh's Beam Model
In the Rayleigh beam model, the effect of shear deformation is neglected in Timoshenko's beam model. Unlike the Timoshenko and Shear beam models, there is only one dependent variable for the Rayleigh beam. Therefore, the potential and kinetic energies can be rewritten as

2.13
Again following the same procedure outline before, the equation of motion in this case is obtained as The general form of any kind of boundary conditions can be listed as Journal of Applied Mathematics 5

Bernoulli-Euler's Beam Model
In this model both effects of shear deformation and rotary inertia are neglected in the Timoshenko beam model. Like the Rayleigh beam model, there is one dependent variable in this model. Based on above assumptions, the potential and kinetic energies in this case are

2.16
Similarly, the equation of motion is and the boundary conditions are as listed in 2.15 .

Homotopy Analysis Method (HAM)
For convenience of the reader, we will first present a brief description of the HAM 8 . Consider the following nonlinear homogeneous differential equations: where N are nonlinear operators, t denotes the independent variable, and W t are unknown functions. Liao constructed the so-called zero-order deformation equation as follows 8 : where q ∈ 0, 1 is an embedding parameter, , h t are nonzero auxiliary parameter and function; respectively, L is an auxiliary linear operator, W 0 t are initial guesses of W t , and φ t; q are unknown functions. As q increases from 0 to 1, the φ t; q varies from the initial approximation to the exact solution. In other words, φ t; 0 W 0 t and φ t; 1 W t . Differentiating once more from 3.2 with respect to the embedding parameter q and then setting q 0, we obtain the first-order deformation equation as which gives the first-order approximation for the W t 8 . The higher-order approximation of the solution can be obtained by calculating the mth-order m > 1 deformation equation where W m t is called the mth-order derivative of unknown function φ t; q .

Timoshenko's Beam Model
For convenience, the following nondimensional variables are used: In order to obtain the nonlinear natural frequency, we start to apply the Galerkin method on which the deflection for the beam with simply supported boundary conditions is expressed as As a result 2.10 can be written as follows: The beam centroid is subjected to the following initial conditions: Note that the free vibration of a system without damping is expressed by a kind of periodic function which has to satisfy the initial conditions. For our case this function simply is cos nω NL t , n 1, 2, 3, . . ..

3.10
It is well known that for a Timoshenko's beam with simple ends there are two natural frequencies. In other words, there are two sinusoidal modes of different frequencies corresponding to the same spatial mode. This is the most important difference between Timoshenko's model and the other beam models with simple end conditions. This phenomenon is also reported by many researchers 9-11 . Thus, the initial guess for W t is as follows: To satisfy the initial conditions given in 3.9 , we apply those conditions into 3.11 which yields the following expressions for the A 1 and A 2 as 3.12 Now we use 3.11 in conjunction with HAM to get the natural frequencies and the beam deflection. To construct the homotopy function, the auxiliary linear operator is selected as

3.13
For convenience and simplify, from 3.7 , the nonlinear operator is defined as

3.15
In our case, to obtain the first-order approximation, 3.3 for the function of φ t; q can be expressed as 3.17 where B's are given in the appendix, equation A.1 . These solutions should obey the general form of the base function. This means that the coefficients of the secular terms cos ω NL1 t , cos ω NL2 t must be zero, that is, B 7 0 and B 8 0. This provides two algebric equations which yield the nonlinear frequencies: 3.18 and solving 3.17 , the W 1 t is obtained as follows: where A and B are obtained by applying initial conditions. Coefficients A, B and D's have been presented in the appendix, equation A.2 . Thus, the first-order approximation of the W t becomes as

Shear, Rayleigh's, and Bernoulli-Euler's Beam Models
The second-order differential equations of motion for either of shear S , Rayleigh's R , or Bernoulli-Euler's B beams using Galerkin's method are the same and is

Postbuckling Load-Deflection Relation and Critical Buckling Load
In order to obtain the postbuckling load-deflection relation one can set all time derivative terms in 3.7 and 3.21 equal to zero which yields the following. Rayleigh's and Bernoulli-Euler's beam models, N π 2 4η 2 4 W 2 .

4.2
Neglecting the contribution of W in the above equations, the critical buckling load can be determined as follows: Rayleigh and Bernoulli-Euler Beam Models:

4.3
Shear and Timoshenko's beam models: N cr π 2 η 2 λπ 2 . 4.4 As it is seen the postbuckling and critical buckling load predicted by Timoshenko model is the same as the one from shear model. The same trend exists between the Rayleigh and the Bernoulli-Euler beam models.

Results and Discussions
In all computational cases, the shear coefficient is taken as k 5/6 and the axial force N 0 unless mentioned otherwise. In order to demonstrate the effectiveness of the HAM, the natural frequencies of a simply supported beam are compared with other existing results in literatures. To do this the nondimensional nonlinear frequencies are calculated, and the results are listed for different η and a in the Table 1.
It can be observed that there is an excellent agreement between the results obtained from HAM and those reported by Rao et al. 6 . It should be noted that only the first natural frequency of a Timoshenko's beam is reported in 6 and in a rather numerical way using FE method. It is clear that Bernoulli-Euler's beam predicts the nonlinear natural frequency higher than the one obtained by other theories. Also the nonlinear natural frequency obtained by Timoshenko's beam theory is lower than the one obtained by other theories. The convergence and accuracy of the presented analytical solution is validated against variation of nondimensional amplitude ratio of the beam center versus t using the Runge-Kutta method. Figure 2 illustrates comparison between these results. They are found to be in a very good agreement. Now, after being satisfied in obtaining results given by this method, in the next step we try to investigate the effect of different parameters on the natural frequencies of four abovementioned beam models.
The variation of the nonlinear frequency versus nondimensional amplitude for all different models and for two different slenderness ratios η using HAM is shown in Figure 3. For thin beams and based on the three employed beam models, namely, Bernoulli-Euler's, Rayleigh's and shear beams, almost no difference is seen for different nondimensional amplitude. However, due to the nature of Timoshenko's simply supported beam, there are two ω NL variations. In the region 0 < a < 5 as it is seen from Figure 3 a , the ω NL1 represents a linear fashion variation whereas the ω NL2 variation remains almost unchanged. In the region a > 20 as it is seen from the same figure, variation of ω NL1 follows a rather constant trend and ω NL2 follows a linear fashion variation. For the intermediate region, that is, 5 < a < 20 frequencies of both ω NL1 and ω NL2 keep changing in a rather nonlinear fashion. As can be seen from Figure 3 a , for thick beams this interval occurs in small amplitude. In other words, for thin and thick beams though in each amplitude both nonlinear frequencies ω NL1 , ω NL2 are involved in beam's response; however, one of them is mostly affected by variation of the nondimensional amplitude. Figure 4 shows nondimensional natural frequency ω NL1 for Timoshenko's beam model versus small nondimensional amplitude for different slenderness ratio and different beam models. It can be seen that the natural nonlinear frequency obtained from four different beam models yields the same value provided that the beam's slenderness ratio at least 100. This finding is also reported by 11 in which only the linear analysis was conducted.
The variations of the nonlinear frequency ratios versus nondimensional amplitude of Timoshenko beam for different slenderness ratios are depicted in Figure 5. As it is seen from these figures the point of shift from flat to nonflat and vise versa of ω NL increases by increasing the beam slenderness ratio. Variation of the nondimensional postbuckling load N versus nondimensional amplitude is shown in Figure 6. As it is seen the postbuckling load predicted by all four different models are the same for rather a slender beam η 50 . On the other hand for the thick beams η 10 , the differences between Timoshenko and/or Shear beam models with Rayleigh's and/or Bernoulli-Euler's beam models are more pronounced. Furthermore, for both cases the postbuckling load increases as amplitude a increases.
Variation of the critical buckling load N cr versus slenderness ratio for the problem under investigation is shown in Figure 7. It can be observed that for η > 12, that is, thin beams, all four models give the same answer for the critical buckling load.

Conclusion
The natural frequencies, deflection, postbuckling, and critical buckling loads of four approximate models of Bernoulli-Euler's, Rayleigh's, shear, and Timoshenko's beams having nonlinear behavior and simply supported end conditions have been discussed in the present survey in detail. A first-order approximation of HAM is used primarily to get the analytic expressions for the nonlinear frequencies and consequently the related deflection of neutral axis of such beam. Based on derived results, the followings are concluded.
1 For a wide range of cases having large amplitude vibration, HAM seems to be very efficient in handling such problems.
2 For thin and thick beams with simple end conditions, there are two nonlinear natural frequencies. Depending on the value of amplitude, only one of these two nonlinear natural frequencies is mostly affected by variation of the amplitude. There is a primary interval of amplitude in which the lower frequency varies in a linear fashion, and the higher frequency remains almost constant. Contrary to this trend, at the ending interval the frequency variations change their character; that is, the lower frequency is almost constant and the higher frequency varies in a linear fashion. Moreover, in the intermediate interval both frequencies increase in a nonlinear fashion. It is also seen that for the thick beams this intermediate range shifts to the left; that is, the primary interval shrinks further. In addition, it is shown that the point of shift from flat to nonflat part of ω NL curves and vise versa increases by increasing the beam slenderness ratio.
3 The Bernoulli-Euler beam model predicts the nonlinear natural frequency higher than the other theories. Whereas, the nonlinear natural frequency obtained by Timoshenko's beam theory is lower than the other theories.
4 It can be seen that the natural nonlinear frequency obtained from four different beam models yields the same value for the beams with slenderness ratio of at least 100.
5 The shear effect is always more dominant than the rotary effect for a given geometry and material in predicting postbuckling and critical buckling loads. Therefore, it is highly recommended to use either of the shear or the Timoshenko models for a thick beam. Nonetheless, in the postbuckling and critical buckling loads analysis