Higher-order nonlinear vibration analysis of Timoshenko beams by the spline-based differential quadrature method

Higher-order nonlinear vibrations of Timoshenko beams with immovable ends are studied. The nonlinear effects of axial deformation, bending curvature and transverse shear strains are considered. The nonlinear governing differential equations are solved using a spline-based differential quadrature method (SDQM), which is constructed based on quartic B-splines. Ratios of the nonlinear to the linear frequencies are extracted and their variations with the ratio of amplitude to radius of gyration are examined. In contrast to the well-recognized finding for the nonlinear fundamental frequency of beams, some higher-order nonlinear frequencies decrease with the increase of ratio of amplitude to radius of gyration.


Introduction
Nonlinear vibrations of beams have been of interest to researchers for many years since they are found in applications of various structural components, from micro-systems to members in conventional mechanical and structural engineering.The majority of the available publications [1][2][3][4][5][6][7][8][9] are based on the classic Bernoulli-Euler beam model which has been deemed capable of evaluating the low order modes and frequencies of slender beams.In recognition of the transverse shear deformation and rotary inertia for stubby beams, Timoshenko beam theory has been employed by some researchers in the nonlinear vibration analysis of beams.
It seems that Rao et al. [10] first examined the nonlinear vibration of Timoshenko beams by setting up a finite element formulation.In their analysis, an improper linearization of the nonlinear axial strain resulted in underestimated nonlinear fundamental frequency.Sarma and Varadan [11] developed a Ritz finite element formulation to study the nonlinear fundamental frequency of a hinged-hinged beam.They reported the increase of the nonlinear fundamental frequency considering the contribution of second order nonlinear curvature.In the case of slender beams, their results were larger than the available exact solution.Furthermore, the use of reduced integration therein during the construction of nonlinear equations for alleviation of shear-lock phenomenon further exacerbated the "hardening" of the nonlinear fundamental frequency.Lin and Tsai [12] also developed a finite element formulation taking into account the nonlinear effects of bending curvature and shear strain in addition to the nonlinear axial strain.In their analysis, the nonlinear effects other than the nonlinear axial strain led to the decrease of the nonlinear fundamental frequency.Adopting a space average approach, Foda [13] studied the nonlinear fundamental frequency of simply supported Timoshenko beams with multiscale method.Similar results were obtained, i.e., the inclusion of transverse shear deformation and rotary inertia brought about decrease of the nonlinear fundamental frequency of Timoshenko beams.More recently, Zhong and Guo [14] established the nonlinear dynamic equations of Timoshenko beams and solved for the nonlinear fundamental frequency using the differential quadrature method.In their analysis, the nonlinear effects of bending curvature and of transverse shear deformation as well as the nonlinear axial deformation were taken into account.They also found that the nonlinear effects of nonlinear transverse shear deformation and bending curvature gave rise to the "softening" of the nonlinear fundamental frequency.An examination of the cause indicates that the conflicting observations are ascribed to the different kinematic relations used in the respective analysis.
The use of Timoshenko beam theory allows the investigation to higher-order modes and frequencies of beams.However, to the best knowledge of the authors, no attempt has so far been made for higher-order vibration analysis of Timoshenko beams.Even for Bernoulli-Euler beams, very limited study has been carried out in regard to the higher-order nonlinear modes and frequencies.The only available work was from Benamar et al. [15].With the harmonic-balance method, they found the first three nonlinear mode shapes of simply supported and clampedclamped Bernoulli-Euler beams.They also reported that the higher-order mode contributions and frequencies increase at large amplitudes.
The differential quadrature method (DQM) [16,17] has been used to perform structural analysis for nearly twenty years.It excels at solution of differential equations over relatively regular geometric domains.Although the employment of Chebyshev-Gauss-Lobatto grid points enables the DQM to cope with various complex problems, the limitation of the number of grid points in the discrete domain still exists.Experience indicates that it is impossible to achieve convergent solution of higher-order nonlinear modes of Timoshenko beams by the conventional differential quadrature method.In contrast, the newly developed spline-based differential quadrature method [18,19], without restriction on the number of grid points, has exhibited edges in the solution of nonlinear initial value problems [20] as well as nonlinear boundary value problems [21].Experience indicates that B-splines of two orders higher than that of differential equations are more suitable for construction of the differential quadrature method for boundary-value problems in light of efficiency and accuracy.In this paper, quartic B-splines are used to construct the spline-based differential quadrature method.The higher-order nonlinear vibrations of Timoshenko beams are then studied.

Spline-based differential quadrature
To construct the basis functions with quartic B-splines for differential quadrature, a set of uniformly spaced nodes are selected in a normalized interval [0, 1], where h is the length of every subinterval.The normalized quartic B-spline function [22] is given by where 5 j are the binomial coefficients, and ∆(x) is the unit step function.Apparently it is a piecewise polynomial that covers six consecutive segments only.To construct a global interpolation function over the normalized interval, extra nodes outside the interval [x 0 , x N ] are usually needed to meet boundary condition requirements, i.e.
A typical spline interpolation function over the given interval can be expressed as where Φ j (x) are usually given in terms of a combination of translated and scaled spline functions.The interpolation condition requires that Φ j (x) satisfy the cardinal condition at every node, i.e. x/h Fig. 1.A quartic cardinal spline function Φ 0 (x).
To this end, three auxiliary spline interpolation functions are constructed, where With the local non-zero property of ϕ 4 (x), all the terms but the one containing y j on the right side of Eqs ( 6) to (8) are eliminated.In consequence, the quartic cardinal spline interpolation function is obtained as Hence, It can be shown that Eq. ( 5) is satisfied at every node.Figure 1 illustrates the quartic cardinal spline function Φ 0 (x).Since the extra nodes outside the interval are often cumbersome to handle, non-integral nodes within the interval in the vicinity of the two ends of the interval are introduced instead, i.e.
The function values at the non-integral nodes are given as where where ) are constants also given in the Appendix.In the same manner, the extra nodes outside the right end of the interval are expressed as Now, the quartic spline interpolation function, without extra outside nodes, is rearranged as follows where Thus, the weighting coefficients in the spline-based differential quadrature are given in explicit form as The localized non-zero nature of splines results in banded weighting coefficient matrices for derivatives.Meanwhile, the following symmetric properties exist among the weighting coefficients , N.

Nonlinear vibrations of Timoshenko beams
A beam with uniform cross-section and immovable ends is considered.Following [14], the kinematic relations adopted for the beam are where u, w and θ are the axial displacement, deflection and rotation, respectively; ε, κ and γ are the axial strain, curvature and the cross-section rotation, respectively.The governing differential equations for nonlinear vibrations of Timoshenko beams with immovable ends are derived via the Hamilton principle, which are where L, A, ρ, I, E, G and k represent the length of the beam, cross-section area, mass density per unit volume, moment of inertia, Young's modulus, shear modulus and shear modification coefficient, respectively.The subscripts x and t denote the partial differentiation with respect to the axial coordinate x and the time coordinate t, respectively.The reader may refer to [14] for the derivation details.
For nonlinear vibration analysis of beams, there are usually two options for simplifying the governing equations.One is the separation of space and time variables, which is believed to be valid for simply supported beams [23].The other assumes that a point of maximum amplitude exists during the vibration and it is also the point of reversal of motion for every point of the beam, which is applicable to beams with other boundary conditions in addition to simply supported condition [24][25][26].The latter option is chosen herein to simplify the governing equations of Timoshenko beams.Namely, it is assumed that at the point of motion reversal, there exist where ω is the nonlinear frequency.
Substituting the above expressions into the governing equations and introducing the following dimensionless parameters the dimensionless form of the governing equations can be written as where a is the maximum vibration amplitude of the beam and the radius of gyration is given by Applying the spline-based differential quadrature rule to Eqs (31) to (32) and invoking the boundary conditions, a set of nonlinear algebraic equations are produced.An iterative procedure [27] is employed to extract the eigenfrequencies.

Numerical results
In all computations, the Poisson's ratio is 0.30 and the shear modification coefficient is 5/6.A benchmark problem -the fundamental frequency of a simply supported slender beam is investigated first using the present nonlinear Timoshenko beam model, for there exists the analytical solution of the nonlinear to the linear fundamental frequency ratio on the basis of Bernoulli-Euler beam model [5,28] where ω * 1L is the analytical solution of the linear fundamental frequency of the beam.In the meantime, the convergence of the spline-based differential quadrature is examined.The relative error of the nonlinear fundamental frequency is defined as which is plotted in Fig. 2 against the number of nodes used in the analysis.It is seen that the relative error is kept below 0.2% if the number of nodes is more than 25.Not only the Timoshenko model is verified but also the efficacy of the spline-based differential quadrature is demonstrated.In addition, the stability of the spline-based differential quadrature is also demonstrated in Fig. 2, where the results of fundamental frequency are virtually invariant for large N .
Table 1 gives the first nonlinear mode shapes of simply supported beams with the increment ratio of amplitude to radius of gyration for L/r = 100 and L/r = 20, respectively.It is found that the mode shapes vary insignificantly with the increment of amplitude ratio for the slender beam (L/r = 100), but significantly for the stubby beam (L/r = 20).It is observed that the largest deviation is merely 0.1% for L/r = 100, but 2.5% for L/r = 20 when the ratio of amplitude to radius of gyration increases from 0.1 to 2.0.Thus, for stubby beams (beams with small slenderness ratio), the variation of mode shapes is noticeable with the variation of ratio of amplitude to radius of  gyration even for simply supported beams.As a result, separation of space and time variables is inappropriate for nonlinear vibration analysis of stubby beams, but is applicable to the vibration analysis of simply supported slender beams with moderately large amplitudes.
Tables 2-4 show the ratios of the first six nonlinear frequencies to the corresponding linear frequencies ω NL /ω L with various slenderness ratios L/r and amplitude ratios a/r for simply supported beams, clamped-clamped beams and beams with one end clamped and the other simply supported, respectively.
In Table 2, the first six ratios of nonlinear frequencies to linear frequencies for simply supported beams are presented.It is seen that the ratios are all larger than unity for slender beams.The ratio increases with the ratio of the amplitude to radius of gyration, and the chosen iterative algorithm runs into difficulty in achieving the convergence of higher-order frequencies.The increase of number of nodes makes the algorithm more sensitive to the selection of initial guest of eigen-mode vector.On the other hand, it becomes even difficult to converge for the higher-order frequencies of stubby beams despite the increase of number of nodes.The results in the table are obtained when N is around 80. It is noted that convergent results cannot be obtained for every chosen node number.A trail and error strategy is therefore adopted to gain a convergent solution.In addition, a distinct discovery is that the fourth and the fifth nonlinear frequencies of the stubby beam are less than their linear counterparts and the ratios Ω 4NL /Ω 4L and Ω 5NL /Ω 5L decrease monotonically with the increase of amplitude to radius of gyration.This is certainly attributed to the intricate nonlinear effects of axial strain, bending curvature and shear strain.It proves difficult to identify the direct cause that leads to the decrease of some nonlinear frequencies since some higher-order frequencies of stubby beams are still found larger than their linear counterparts, as in the case of the sixth nonlinear frequency in Table 2. Similar situations arise for clamped beams and the beams with one end clamped and the other simply supported, as shown in Tables 3 and 4. The mode shapes of a stubby beam with one end clamped and the other simply supported are illustrated in Fig. 3.In a nutshell, the vibration amplitude acts as a magnifier that controls either the increase or the decrease of the nonlinear to the linear frequency ratios.

Conclusion
Higher-order vibrations of Timoshenko beams have been studied using a spline-based differential quadrature method which is constructed with quartic B-splines.A quartic cardinal spline function is given and the explicit weighting coefficients for derivatives are obtained.Attempts have been made to find the first six nonlinear frequencies of clamped beams, simply supported beams and beams with one end clamped and the other simply supported.Although the convergence of the higher-order frequencies proves difficult with the increase of vibration amplitude, especially for stubby beams, the variation tendency of the nonlinear frequencies has been identified.In particular, some higher-order frequencies are found lower than their linear counterparts.Further research is needed to probe into the cause of the decrease of some nonlinear frequencies of stubby beams with the increase of ratio of vibration amplitude to radius of gyration.The coefficients in Eqs ( 16) and ( 17

Fig. 3 .
Fig. 3. First six nonlinear mode shapes of a beam with the left end clamped and the right end simply supported (L/r = 20, a/r = 0.3).

Table 1
Nonlinear first mode shape for simply supported beams

Table 2
Ratio of nonlinear frequencies to linear frequencies for simply supported beams

Table 4
Ratio of nonlinear frequencies to linear frequencies for beams with one end clamped and the other simply supported