Choosing the Best Functions for the Rayleigh-Ritz Vibration Analysis of Beams

The Rayleigh-Ritz Method (RRM) is used extensively for the vibration analysis of structures. The accuracy depends on the assumed functions. In this work, we include both sine and cosine functions in the admissible functions, which is the first time to represent symmetric and anti-symmetric modes. Several different groups of functions are examined and compared for the accuracy of the resulting natural frequencies, and for the overall mode shape error norms calculated with respect to the known exact solutions. It is concluded that a set that combines low order polynomials, odd cosine and odd sine functions, or, even cosine and even sine functions, is more likely to yield the best accuracy and convergence of both frequency and mode shapes for a general beam structure.


Introduction
The Rayleigh-Ritz method (RRM) [1]- [4] has been used extensively for structural vibration problems.A complete presentation of the method and its use was published [5].The main factor in this method is the selection of the functions that are used to represent the system.Over the years several types of solution were offered: polynomial functions, trigonometric functions and combinations of these.Monterrubio and Ilanko [6] review these functions and consider a combination of low order polynomials and cosine functions.In this work several different sets of functions are compared for the problem of beam vibrations.The comparisons are presented for the two characteristics of the vibration problem: the natural frequencies (eigenvalues) and vibration modes (eigenvectors).For the first one can calculate a representative relative error, and the second will be assessed by defining a normalized error norm of the difference between the approximate vibration modes and the exact modes.
The study is performed for the classical beam theory of Bernoulli-Euler, for which the exact vibration characteristics are known.The enforcement of boundary condition is performed by adding penalty matrices which account for all the 4 possible end restraints [5]- [8].
In past research the sets of functions that were used for the analysis include: (a) polynomial functions that exhibit fast convergence but are also prone to ill-conditioning, (b) polynomial functions that are improved by performing the Gram-Schmidt orthogonalization procedure [9], (c) Fourier cosine function with additional polynomial terms [10], and Fourier sine series with additional polynomial terms [11].In all these studies the attention was focused on the vibration frequencies with no account for the vibration modes approximations.More recently, Koo et al. used polynomials to study the natural modes of a telescopic boom system with multiple sections.They obtained the first-mode natural frequency with an error of 4.6 % compared to FEM [12].Bao et al. used two sine functions and several cosine functions to conduct a vibration analysis of nanorods [13].They showed the accuracy of the results with eight terms, but the convergence of the results is not monotonic.Babaei investigated longitudinal vibration responses of axially functionally graded optimized MEMS [14].They used sine functions or cosine functions depending on the boundary conditions considered.Ozbasaran evaluated the buckling analysis of the I-section prismatic beam-columns [15].They considered three different series of admissible functions, namely power, trigonometric, and exponential trial functions.They used different forms for different restraint configurations.It was shown that the trigonometric series give better results compared to the power and exponential series forms.Yang et al. presented a general approach for the free vibration analysis of a circular cylindrical shell resting on elastic foundation [16].They used exact beam functions as admissible functions.The results showed good convergence and accuracy with 12 terms but only for the (1,1) mode.Mazanoglu proposed part by part usage of Timoshenko and Euler-Bernoulli beam theories for obtaining natural frequencies of the non-uniform beam [17].The author also presented convergence tests to determine the proper function among the simple admissible shape functions, such as polynomial functions and trigonometric functions.It is shown that using 12 terms gives good accuracy, and trigonometric functions are superior to polynomial functions.Fakhar and Hosseini-Hashemi studied the static bending and free vibration behaviour of Euler nanobeams using the combination of polynomial and trigonometric functions [18].They showed that using more terms gives more accurate results without ill conditions.
Baruh and Tadikonda [19] have introduced the requirement for the admissible functions to satisfy the additional Complementary Boundary Conditions (CBC).These are defined as the relations that arise due to terms in the boundary expressions in the variational form of the equations.It was shown that if the set of admissible functions violate the CBC it is equivalent to imposing additional constraints on the solution that will result in slow convergence, or no convergence to the correct solution at all for some cases.This was taken care of in the formulations in [5] by adding terms that enable the enforcement of the CBC at the end of the beam.Most of the researchers in the past have used the RRM for plate vibration problems, however they introduced one dimensional (beam) functions in two directions to represent the assumed admissible surface of the plate.However, even for beams, unless the CBC are satisfied, sometimes it is not possible to get convergence to the correct solution as the function sets used may be over constrained [20].
In the present work, the emphasis is put on beam vibration problems.We include both sine and cosine functions in the admissible functions to represent symmetric and anti-symmetric modes.It is shown that including both sine and cosine functions improves the accuracy of not only the natural frequencies but also of the mode shapes.
In the next section a brief summary of the method is presented, followed by five function groups that are tested.Then the derivation of the stiffness, mass, and penalty matrices for the analysis is given.Numerical results are given for several modes and boundary conditions representative of all the possible cases.

The Rayleigh-Ritz Method (RRM) for Beams
For the sake of completeness, some derivations for beam vibration problems and the Rayleigh-Ritz method (RRM) are taken from [5] and presented here.
For Bernoulli-Euler beams the equation of motion is where, E:Young's Modulus, I:Second moment of area, ρ: density, A: cross-sectional area, x:the axial coordinate, y:lateral displacement, ω: circular frequency of the beam respectively.We obtain the terms for the total potential energy and kinetic energy as In the RRM the functions () are assumed as the sum of a series product of undetermined coefficients and admissible functions   (ξ), with ξ = / as the normalized coordinate along the member, given by In the RRM the natural frequencies are obtained as where   ,   are the maximum values of the potential energy and the kinetic energy function, respectively.As the Rayleigh Quotient (RQ) given by equation ( 5) gives an upper-bound to the natural frequencies, to obtain the best estimate of the natural frequencies we minimize the RQ w.r.t the undetermined displacement coefficients, G's.
for i=1, 2, .., N.This results in the following system of equations, with [, ] =  2       (9) The enforcement of the three cases of boundary condition, i.e. free (F), simply supported (S) and clamped (C), at the two ends of the beam is done using the penalty formulation (artificial springs) as derived by Ilanko et.al. [5].The translational restraints (p1 and p3) and the rotational restraints (p2 and p4), shown in Figure 1, have a contribution to the system potential energy, and it results in the addition of the following matrix [P] and the values of pi, i=1...4 are taken to satisfy the end restraint, for example, p1 = p2 = 0 for free, p1 = 0 and p2 = ∞ for simply supported, and p1 = p2 = ∞ for clamped conditions respectively.(We can also set a numerical value for a specified spring stiffness for partial restraint).Then Eqn.(7) will have the final form The solution of Eqn.(7) gives the normalized frequencies of vibration for the beam and the coefficients   that can then be used in Eqn.(4) to find the vibration modes.

Function Groups
For the vibration analysis of beams using the Rayleigh-Ritz method the functions used will have to be able to satisfy all the standard boundary conditions that are at the two ends: Free (F), Simply Supported (S), Clamped (C), and Guided (G).The deflected shape of the beam at these locations is required to have zero deflection (S and C), zero slope (C and G), and zero curvature (S, F), or zero shear force (F, G) or a combination of these.The trigonometric functions cosine and sine can be used for some of the combinations and provide good approximations for the mode.Obviously for some boundary conditions they will give exact results, for example sine series for beams with simply supported ends and cosine series for beams with gliding ends.In some cases the use of these functions will make the compliance with the required boundary values more complicated, and a larger number of functions will be required for good results.These considerations were taken in the groups of functions that are proposed in this work.
As the terms in the stiffness and mass matrices ([K] and [M]) involve integrals of multiplications of the selected functions, the use of the cosine and sine functions will have the benefit of orthogonality, such that for some of the integrals we shall have zero values, and in some, the resulting matrices will be diagonal.This will speed up the convergence and reduce the influence of numerical and round-off errors in the calculations.
In order to be able to satisfy the boundary conditions it was shown in Ref. [5] that we need to have the following three functions: These three functions are needed to enable translations and rotations at the two ends of the beam segment.In addition to these three functions it is possible to add more functions that will improve the accuracy of the natural frequencies and the mode shapes of the beam.
The choice of the additional functions is wide and several possibilities are examined in this work and are grouped as follows.All the groups will have the first 3 function φ1(ξ), φ2(ξ), and φ3(ξ) and some additional functions.
(a) Function Group 1 -FG1: In this group we add the Fourier cosine series with N1 terms: This series of functions was used in the calculations in Ref. [5].
(b) Function Group 2 -FG2: In this group we add the Fourier sine series with N2 terms: This group is similar to FG1, however it is expected to yield faster and more accurate results for cases with hinged ends.
(c) Function Group 3 -FG3: In this group we add the Fourier cosine and sine series with N3 terms: This group is combining the previous two.

(d) Function Group 4 -FG4:
In this group we add only the odd terms in Fourier cosine and sine series with N4 terms: Using these functions we still preserve the orthogonality of the trigonometric functions as was the case for the FG1 and FG2 cases.Mass matrix has off-diagonal terms only in the first three rows and columns which correspond to the terms involving the linear and square polynomials [6], and this helps to zero most of the offdiagonal terms.
(e) Function Group 5 -FG5: In this group we add the only the even terms in Fourier cosine and sine series with N5 terms: In this group also, all the trigonometric functions in the group are orthogonal to each other, as was the case for the FG1 and FG2 groups.

The Stiffness, Mass and Penalty Matrices
In this section the matrices [K]j, [M]j and [P]j for the five function groups, j=1...5, (FG1 to FG5), are given as follows.All the matrices are symmetric so only the upper triangle will be given.The matrices are subdivided into 4 submatrices as follows where Nf is the number of functions taken in each function group (FG).The explicit terms in these matrices are given in Appendix 1.

Orthogonality of the Admissible Functions
The terms in the [K]j and [M]j involve integrations and differentiation of the assumed functions, as follows where those for the penalty matrices involve evaluations of the functions at the boundary points (Eqs.11-14).As was shown, in some function groups the functions are orthogonal and this will result in simpler K and M matrices, as follows: (a) FG1 -K matrix is diagonal, MBB is diagonal (b) FG2for this group we will have more complex terms, and some of the values in the penalty function related to the derivatives of the function will not be equal to zero.(c) FG3 -In using both sine and cosine terms (seemingly more complete expansion) the orthogonality of many terms is violated and many new terms will be present in all the matrices.This is expected to cause convergence difficulties and increase the relative errors.(d) FG4 -In this group the sine and cosine functions are taken as the orthogonal terms only, and one will expect better performance as the boundary conditions could be satisfied by fewer terms.(e) FG5 -This group is similar to the fourth group, and may have better performance for certain combinations of boundary conditions.
The accuracy in the computation of the eigenvalues and eigenvectors is affected significantly by the structure of the matrices involved.The better FG which will be recommended is the one that result in smaller relative numerical errors in both the vibration frequencies and modes.

Numerical Results and Discussion
The natural frequencies and mode shapes are obtained using MATLAB.The comparison of the numerical accuracy of the various groups will be done using two comparison criteria: 1. Relative error of the natural frequencies, compared to the exact solutions available in the existing literature [5], [21], and 2. The magnitude of the error norm in the mode shape taken as where   () is the high accuracy mode shape for the particular set of boundary conditions given by Blevins [21], and N is the number of functions used for the solution.Both functions are normalized such that the maximum value of the deflection in a particular mode is equal one, so that one can obtain meaningful results for the norm.In the results that are presented, for each frequency the error norm is calculated with respect to the same exact known mode shape, and thus can be used for comparison within that group.In the numerical results for the different groups the number of functions that were used was the same, i.e.N = Ni + 3, with N1 = N2 = N4 = N5 = 2*Nf, and N3 = Nf.Using the penalty function to enforce the restraints results in the need for minimal number of functions in the analysis in order to obtain reliable results as some are used to satisfy the boundary conditions.For each of the cases presented below there is a different number, as the number of restraints is different.So, in the following presentation of the results this minimum number of functions is given for every frequency.The tables include the meaningful results for the cases and are given until the relative error is very small.For several of the cases we can identify symmetric and anti-symmetric modes and these are also pointed out in the tables.The type of functions in each column in the tables, is given where C stands for cosine functions, and S stands for sine functions.The letters O and E indicate odd functions only, and even functions only, respectively.
(a) Free -Free Beam Tables 1 and 2 give the results for the F-F case where no penalty matrices are needed.It can be seen that the better convergence among the five FG, for both the value of the frequency and the error norm of the modes, is not the same for all the modes.It alternates between the FG1, FG4, and FG5 combinations, for the frequencies and between the FG1, FG2, and FG5 combinations for the normalized modal error, dependent also on the mode being symmetric or antisymmetric for the particular frequency.

(b) Clamped -Clamped Beam
The results for this case are given in Tables 3 and 4. For this case results from the functions in FG4 are the best for both frequency and modal normalized errors.The results from FG1 and FG2 groups are equivalent to the best results from FG4, for the symmetric (FG2) and anti-symmetric (FG1).

(c) Clamped -Simply Supported Beam
The results for this case are given in Tables 5 and 6.This case in not symmetric in terms of boundary conditions as were the first two.Then, we expect that since the modes cannot be cast into symmetric and anti-symmetric classes, there is an advantage to the function groups that will include both sine and cosine functions (FG3, FG4, and FG5).In the results we can observe the FG3 performs the best for the lower frequencies and modes, but for the higher modes FG4 has the advantage.This is explained by the fact that not all the functions in FG3 are orthogonal to each other, and as was the case for other studies, in such cases numerical accuracy is affected significantly from the large number of off diagonal terms in the matrices, when the number of function in the group is increased.The results for this case are given in Tables 7 and 8.This case is also not symmetric and the results obtained from FG3 are the best, subject to the restriction in the number of functions.Overall, for this case FG4 yields the best performance.
As was shown in the numerical results, in each case the advantage is to either the FG4 or FG5 group.For Free-Free beam, FG1 also give as good results as FG4 or 5 does for the natural frequencies.However, for mode shapes, FG5 always gives the best results.The matrices that are used for both FG4 and FG5 have small differences, and one can choose which group yields better results by performing the calculation for both groups using a small number of members in each group (say 2 sine and 2 cosine terms).Comparing the results of these two analyses would reveal the better choice, and then the computations can be continued using that group.For plate problems, this strategy may yield much faster convergence and higher accuracy, and will be explored in future work.

Summary and Conclusions
In the vibration analysis of beams using the RRM the accuracy of the solution is highly dependent on the admissible functions that are used in the computations.Two factors are important in this regard: the ability of the set to satisfy the boundary conditions, and the number of functions that are required to obtain high accuracy.It is shown in this work that choosing a group of functions that are capable to model all the possible combinations will benefit the most.This is achieved by using the combined odd sine and odd cosine functions (FG4), or the combined even sine and even cosine functions (FG5), which is the first time both functions are included simultaneously in the admissible functions.In both, all the functions in the group are orthogonal, so one also preserves the high numerical accuracy that is obtained in the computation that is performed with diagonal matrices.
In this work a new error norm for the mode shapes was also suggested at the first time.The proposed admissible functions give good accuracy also for the mode shapes.It was shown that it has great importance in the evaluation of functions which are used in the RRM, and will also be expanded for the analysis of plates in future work.Appendix: Matrices for Function Groups FG1-5 For all the matrices, the sub-matrices AA that include the first 3 rows and first 3 column of the stiffness, mass and penalty matrices are identical for all the function groups and are given by:

Function group 5 (FG5):
The matrices for this group are the same as for the group 4, but only even values are taken.The P matrices are the same as for FG3 with: () = 2, 4, 6… and ,  = 2,4,6 …

Table 2 .
Error norm of the normalized vibration modes -FF Beam

Table 3 .
Relative error in frequency -CC Bean

Table 1 .
Relative error in frequency -FF Beam

Table 4 .
Error norm of the normalized vibration modes -CC Beam

Table 5 .
Relative error in frequency -CS Beam

Table 6 .
-Error norm of the normalized vibration modes -CS Beam

Table 7 .
Relative error in frequency -CF Beam

Table 8 .
Error norm of the normalized vibration modes -CF Beam