A spectral finite element model for vibration analysis of a beam based on general higher-order theory

The spectral element matrix is derived for a straight and uniform beam element having an arbitrary cross-section. The general higher-order beam theory is used, which accurately accounts for the transverse shear deformation out of the cross-sectional plane and antielastic-type deformation within the cross-sectional plane. Two coupled equations of motion are derived by use of Hamilton’s principle along with the full three-dimensional constitutive relations. The theoretical expressions of the spectral element matrix are formulated from the exact solutions of the coupled governing equations. The developed spectral element matrix is directly applied to calculate the exact natural frequencies and mode shapes of the illustrative examples. Numerical results of the thick isotropic beams with rectangular and elliptical cross-sections are presented for a wide variety of cross-section aspect ratios.


Introduction
Beam members are basic structural components widely used in many types of structures within the fields of mechanical, aeronautical, automobile and civil engineering.Because of the practical importance, a large number of investigators have made efforts to deal with the static and dynamic analyses of the beams.The analyses of slender beams have been performed mostly using the elementary Bernoulli-Euler theory.This theory assumes that the plane sections of the cross-section remain plane and particular to the beam axis.
A more refined beam theory is Timoshenko beam theory which accounts for the shear deformation and rotary inertia that are neglected in the elementary theory.In Timoshenko beam theory the distribution of the transverse shear strain is constant through the beam thickness and therefore requires the shear correction factor to correct the strain energy of deformation.Timoshenko [1] first developed a theory which allowed one to study the vibrational behavior of the thick beams by approximately accounting for transverse shear deformation and rotary inertia, where the cross-section remained undistorted and the in-plane stresses were zero.Mindlin and Deresiewicz [2] calculated the shear correction factor for a variety of cross-sections of the beams.The dynamic shear strain distribution may differ significantly from the parabolic form of the static shear strain distribution.Cowper [3] presented a method to determine the shear correction factor based on the static three-dimensional elasticity theory along with Saint-Venant's static flexure warping function.Hutchinson [4] investigated the behavior of the shear coefficient for a rectangular Timoshenko beam by comparison of the Timoshenko beam solution with a three-dimensional solution for a simply-supported beam.It was found that a shear coefficient would have to be a function of the wave length as well as the aspect ratio.
The limitations of the elementary Bernoulli-Euler beam theory and first-order Timoshenko shear deformation beam theory force the development of higher-order shear deformation beam theories.Stephen and Levinson [5] developed a second-order model for studying the beam vibrations with two governing differential equations having two independent coefficients.The first coefficient was similar to the Timoshenko's shear correction factor, which depended upon Saint-Venant's static flexure warping function, while the second coefficient depended upon the transverse shear stress state.Levinson [6] developed the higher-order shear deformable theory for beams with thin rectangular cross-sections that correctly accounted for the stress free boundary conditions on the upper and lower surfaces as well as the parabolic shear strain distribution through the thickness without the need for a shear correction factor.This was accomplished by expanding the axial displacement to include a cubic distribution through the thickness.This additional displacement was identical to the Saint-Venant's static flexure warping function for thin rectangular cross-sections, but for general cross-section shapes the correct expansion should be an infinite series of transcendental functions.Bickford [7] considered a shear deformation theory for rectangular beams that correctly accounted for the stress free boundary conditions on the upper and lower surfaces of the beam while retaining the parabolic distribution of the transverse shear strain.Murty [8] formulated a third-order beam theory including the transverse shear strain and nonclassical axial stress.In this theory the parabolic transverse shear stress distribution across the thickness of the beam could be obtained using the constitutive relations.Levinson [9] obtained the higher-order beam theory with the fourth-order differential equations, satisfying two boundary conditions at each end of the beam.No shear correction factors were required since the theory satisfied the shear stress free surface conditions on the top and bottom of the beam.Rychter [10] verified the consistency and accuracy of Levinson's theory by imbedding it in the two-dimensional linear theory of elasticity.It was shown that the Levinson's variationally inconsistent theory was more accurate than the variationally consistent theory of Bickford.Heyliger and Reddy [11] used a higher-order theory and noted the inadequacy of the shear correction factor in predicting the natural frequencies.Petrolito [12] presented a finite element solution of uniformly loaded thick beams using the Bickford's theory.Kathnelson [13] developed an improved engineering theory for elastic isotropic uniform beam using the asymptotic analysis of the variational equation of Hellinger-Reissner principle for the three-dimensional theory of elasticity.Wang et al. [14] discussed the derivation and comparison of the Bernoulli-Euler, Timoshenko, and Reddy-Bickford beam theories extensively.Reddy et al. [15] gave solutions to the Levinson beams in terms of the classical beam solutions.Ghugal and Shimpi [16] presented a review of the displacement and stress based refined theories for isotropic and anisotropic laminated beams.Aalami [17] studied the vibrational behavior of a general prismatic beam with an arbitrary cross-section by using a Rayleigh-Ritz energy approach to the general three-dimensional problem.His numerical and graphical results clearly illustrated the presence and importance of both the out-of-plane shear deformation and in-plane deformation for extremely high vibrational modes.Ie and Kosmatka [18] developed a static theory for a general cylindrical or prismatic beam that incorporated both out-ofplane shear and in-plane deformation functions, where these functions were assumed known.In fact, these functions could be determined exactly for simple cross-sections by solving Saint-Venant's bending and flexure problems, and approximately for an arbitrary cross-section by applying either a two-dimensional finite element approach or a power series approach.Kosmatka [19] presented a new beam theory for studying the static and vibration behaviors of the cylindrical or prismatic beam-type structures having an arbitrary cross-section.This general higher-order theory accurately accounted for the transverse shear deformation out of the cross-sectional plane and antielastic-type deformation within the cross-sectional plane.
A review of the literature indicates that the refined theories for dynamic analysis of the beams with non-rectangular cross-sections are rare.Most of the research work related to the dynamic analysis of thick isotropic beams assumes that the in-plane cross-sectional stresses are negligible and the cross-section does not deform in its own plane and comparatively less work is available on incorporating both out-of-plane shear and in-plane deformation in a unified manner.The objective of this paper is to develop an exact spectral element matrix for a general higher-order beam element with an arbitrary cross-section, which has not still been reported in literature.Although several researchers have investigated and elaborated the spectral element method [20,21], it should be noted that there are differences in the way they compute the spectral elements from the way used in this paper.The explicit analytical expressions for the entries of the spectral element matrix are derived in the paper.The spectral element method is often referred to as an exact method as it is based on the exact shape functions obtained from the exact solutions of the element differential equations.The usefulness of the method becomes apparent when higher frequencies and better accuracies of the results are required.In this paper, the spectral finite element formulation for a thick isotropic beam having an arbitrary cross-section is presented.The general higher-order beam theory is employed [19], which considers both out-of-plane shear dependent warping and in-plane deformations.The theoretical expressions of the spectral element matrix for the general higher-order beam element are derived from the exact solutions of the governing differential equations.The application of the derived spectral element matrix to calculate the natural frequencies and mode shapes of the particular beams uses the automated Muller root search algorithm [22].Numerical results are presented and compared with the available solutions in literature to validate the accuracy and efficiency of the current approach.

Mathematical formulation
A prismatic isotropic beam with a general homogeneous cross-section, as shown in Fig. 1, is considered in this paper.A Cartesian co-ordinate system (x, y, z) is defined on the beam, where x and y are coincident with the principal axes of the beam cross-section and z is coincident with the centroidal axis.The current development is further restricted to the study of vibrational behavior in the y − z plane only, and the displacement relations are defined as [19] u 1 (x, y, z, t) = ∂θ(z, t)/∂zψ x (x, y) (1a) where (u 1 , u 2 , u 3 ) are the displacements of a point along the (x, y, z) coordinates, v(z, t) represents the timedependent y direction displacement and θ(z, t) the time-dependent rotation about the x axis.The three functions ψ x (x, y), ψ y (x, y), and ψ z (x, y), which are assumed to be known, can be determined exactly for simple crosssections by solving Saint-Venant's bending and flexure problems, or approximately for arbitrary cross-sections using either a two-dimensional finite element approach or a power series approach [19].
The time-dependent strain-displacement relationships of the beam are defined as ε yy = ∂θ/∂z∂ψ y /∂y (2b) where ε xx , ε yy and ε zz are the normal strains in the x, y and z directions, respectively, and γ yz , γ xz and γ xy are the shear strains.
The general constitutive relations are given as Here, λ and G are defined as Lame's constants.
where E is Young's modulus and ν is Poisson ratio.The total strain energy U of the beam shown in Fig. 1 is given by Substituting σ xx , σ yy , σ zz and τ yz , τ xz , τ xy from Eqs (3a)-(3f) into Eq.(4a) yields The total kinetic energy T of the beam is defined by where ρ is the mass density per unit volume of the beam and the superscript dot represents a derivative with respect to time t.Substituting u 1 , u 2 and u 3 from Eq. (1) into Eq.(5a) results in where the superscript prime refers to a partial derivative with respect to z co-ordinate.The governing equations of motion and boundary conditions for the displacement relations given by Eq. ( 1) can be obtained conveniently by means of Hamilton's principle, which can be stated in the form Herein δU and δT are the variations of the strain and kinetic energies, respectively.Substituting Eq. (2) into Eq.(4b), then substituting the resulting equation and Eq.(5b) into Eq.( 6) and performing the variational operations yields the following governing differential equations of motion for the beam based on general higher-order theory The associated boundary conditions at the ends of beam (z = 0, L) are obtained as follows

Spectral finite element
It can be shown that Eq. ( 7) have solutions that are separable in time and space, and that the time dependence is harmonic.Letting where ω is the angular frequency, V (z) and Θ(z) are the amplitudes of the sinusoidally varying transverse displacement and normal rotation, respectively.Introducing Eq. ( 9) into Eq.( 7), the following differential eigenvalue equations are obtained The solutions to Eq. ( 10) can be expressed as Substituting Eq. ( 11) into Eq.( 10), the equivalent algebraic eigenvalue equations are obtained and the equations have nontrivial solutions when the determinant of the coefficient matrix of Ã and B vanishes.Setting the determinant equal to zero yields the characteristics equation, which is an eighth-order polynomial equation in κ where Therefore the fourth-order polynomial equation for the roots χ = κ 2 must be solved.The solutions can be found in close-form as follows.
Equation ( 12) can be rewritten as where The fourth-order Eq. ( 13) can be factorized as where a 2 1 − 4a 2 + 4λ 1 and λ 1 is a real root of the following cubic equation Then the four roots of the fourth-order Eq. ( 13) can be written as If the discriminant of cubic Eq. ( 14), D, is negative, then Eq. ( 14) has usually three real roots, which can be written as Sign convention for positive shear force Q(z), bending moment M (z), generalized moment M (z) and higher-order moment P (z).
Whichever of the three roots is chosen, the same solutions are obtained.
The general solutions to Eq. ( 10) are given by where In the solution of Eq. ( 13), if any of the χ j 's are zero or are repeated, the solutions to the differential Eq. ( 10) will be modified according to the well-known methods for ordinary differential equations with constant coefficients, for those particular values of χ j .
From Eq. ( 16), only eight of the sixteen constants are independent.The relationship among the constants is given by where According to the sign convention shown in Fig. 2, the expressions of shear force Q(z), bending moment M (z) , generalized moment M (z) and higher-order moment P (z) can be obtained from Eqs ( 8) and ( 16) as follows Θ Fig. 3. Boundary conditions for displacements and forces of beam element.
Refer to Fig. 3, the boundary conditions for displacements and forces of the beam element are, respectively, Substituting Eq. (18a) into Eq.( 16), the nodal displacements defined in Fig. 3

can be expressed in terms of
where {D} is the nodal degree-of-freedom vector defined by e κ2L e κ3L e κ4L e −κ1L e −κ2L e −κ3L e −κ4L t 1 e κ1L t 2 e κ2L t 3 e κ3L t 4 e κ4L −t 1 e −κ1L −t 2 e −κ2L −t 3 e −κ3L −t 4 e −κ4L κ 1 e κ1L κ 2 e κ2L κ 3 e κ3L κ 4 e κ4L −κ Substituting Eq. (18b) into Eq.( 17), the nodal forces defined in Fig. 3 also can be expressed in terms of A i as where {F } is the nodal force vector defined by 19) and (20) gives the following relationship between the nodal forces and nodal displacements where [K] is the spectral element matrix.It should be mentioned that the explicit analytical expressions for the terms of the spectral element matrix could be derived using the symbolic manipulator software such as Maple [23], although the expressions are too lengthy to list in the paper.

Muller root search method
Once the spectral element matrix is obtained, the appropriate boundary conditions for the particular beam under consideration are applied to obtain the frequency characteristic determinant, which equals zero when an eigenvalue of the beam is found.Obviously the resulting frequency characteristic equation is a transcendental function of frequency, which allows an infinite number of natural frequencies to be accounted for.Although more advanced algorithms and more complex procedures may be used to determine the frequency values of the characteristic equation, a simple automated Muller root search method [22] is adopted in the present study to determine all the natural frequencies in a given frequency band.The mode shapes corresponding to the natural frequencies can be found in the usual way by making an arbitrary assumption about one unknown variable of the beam and then calculating the remaining variables in terms of the arbitrarily chosen one.
Let us denote the frequency characteristic determinant by f (ω) and consider the solution of f (ω) = 0 by Muller's method.Muller's method uses a quadratic approximation to the function f by interpolating a quadratic polynomial through the last three computed points and then determining where this curve crosses the ω axis.The following algorithm is an implementation of Muller's method.This algorithm terminates when a prescribed error tolerance on |f (ω)| has been achieved or a maximum number of iterations has been performed.
(1) Pick three distinct values of ω, i.e., ω 0 , ω 1 , ω 2 .Set i = 1 and compute  (3 where ω i+2 is the new approximation of the root ω i+1 .The process is continued until a desired accuracy is obtained.The rate of convergence is high and may be increased if an estimate of the eigenvalue is available.Previously roots may be removed from f (ω) by dividing it by

Numerical results and discussion
Two example problems are considered to demonstrate the accuracy of the spectral finite element model based on general higher-order beam theory.
The first example is of a rectangular cross-section beam shown in Fig. First, the vibration behavior of the Simply-Supported beam is studied for three cross-section aspect ratios b/a.To validate the correctness and accuracy of the proposed formulation, the present results are given in Table 1 and compared with the previously published solutions [19] whenever possible.The numerical results based on third-order shear deformation and Timoshenko beam theories are also shown in Table 1, the shear correction factor is taken as 5/6 for Timoshenko beam theory.It can be seen from Table 1 that the present results are excellent agreement with the existing solutions.As a further validation of the current spectral element model, Aalami [17] predicted the fundamental frequency for a beam with L = 10 m having a square cross-section with a = b = 5 m using a Rayleigh-Ritz energy approach on the full three-dimensional problem, which is 122.09Hz.This value is very close to the present solution 122.62 Hz.
Then, the spectral finite element model of the beam based on the general higher-order theory is used to compute the natural frequencies of the beam with a variety of aspect ratios.Five different aspect ratios, i.e., b/a = 1000, b/a = 10, b/a = 2, b/a = 1 and b/a = 0.5 are considered.In order to investigate the effect of the aspect ratio on the vibration frequencies, the first six natural frequencies of the beam are calculated and the percentage variations are plotted against the aspect ratio in Fig. 5.The percentage variation is defined as percentage variation = 100(ω − ω 0 )/ω 0 ,  where ω 0 is the natural frequency of the beam with the aspect ratio b/a = 0.5.The numerical values of the first six natural frequencies are also shown in Tables 2 and 3.For the Simply-Supported boundary condition, the calculations of the analytical solutions using the expressions given by Kosmatka [19] are also performed.However, the analytical solutions are not displayed in Table 2 because they are complete agreement with the results predicted by the present formulation.It can be noticed from Fig. 5 that the effect of aspect ratio b/a on the natural frequencies of the beam is rather complex, it depends upon the particular boundary condition and mode number.There is a tendency that the aspect ratio has a pronounced effect on the higher frequencies.For the Clamped-Free end condition, increasing the aspect ratio from 0.5 to 1000 yields an increase in the frequency approximately −10.45%, 6.44%, −0.83%, 38.69%, 12.46% and 30.87% for the first, second, third, fourth, fifth and sixth modes, respectively.For the Clamped-Free boundary condition, the first six normal mode shapes of the beam with cross-section aspect ratio b/a = 1000 are computed and the numerical results are plotted in Fig. 6.
The second example is a beam with elliptical cross-section, as shown in Fig. 4. The numerical data for this example is the same as the previous one.The first six natural frequencies of the beam with five different cross-section aspect ratios, i.e., b/a = 1000, b/a = 10, b/a = 2, b/a = 1 and b/a = 0.5 are also calculated and the numerical results are given in Tables 4 and 5.The percentage variations of the natural frequencies are plotted against aspect

2 Fig. 4 .
Fig.4.A rectangular section used in example 1 and an elliptical section used in example 2.

Fig. 5 .
Fig. 5. Percentage variations of first six frequencies versus aspect ratio for (a) Clamped-Free and (b) Simply-Supported rectangular section beams.

Table 1
Comparison of natural frequencies (in Hz) of Simply-Supported rectangular section beam

Table 2
Natural frequencies (in Hz) of Simply-Supported rectangular cross-section beam