Power secant method applied to natural frequency extraction of Timoshenko beam structures

Power secant method applied to natural frequency extraction of Timoshenko beam structures

C.A.N. Dias^*

Group of Solid Mechanics and Structural Impact, Department of Mechatronics and Mechanical System Engineering, University of São Paulo, São Paulo – 05508-900 – Brazil

ABSTRACT

This work deals with an improved plane frame formulation whose exact dynamic stiffness matrix (DSM) presents, uniquely, null determinant for the natural frequencies. In comparison with the classical DSM, the formulation herein presented has some major advantages: local mode shapes are preserved in the formulation so that, for any positive frequency, the DSM will never be ill-conditioned; in the absence of poles, it is possible to employ the secant method in order to have a more computationally efficient eigenvalue extraction procedure.

Applying the procedure to the more general case of Timoshenko beams, we introduce a new technique, named "power deflation", that makes the secant method suitable for the transcendental nonlinear eigenvalue problems based on the improved DSM.

In order to avoid overflow occurrences that can hinder the secant method iterations, limiting frequencies are formulated, with scaling also applied to the eigenvalue problem. Comparisons with results available in the literature demonstrate the strength of the proposed method. Computational efficiency is compared with solutions obtained both by FEM and by the Wittrick-Williams algorithm.

Keywords: exact modal analysis, dynamic stiffness matrix, secant method, power deflation.

1 INTRODUCTION

Using an improved DSM, exact mode shape and natural frequency calculation of Euler-Bernoulli skeleton structures is discussed at length by Dias and Alves [4]. According to the Euler-Bernoulli theory, the beam transverse vibrations are only due to bending deflection. Hence, shear deflection is not taken into account, which restricts the work in [4] to the slender members and low order modes.

For short beams, the Euler-Bernoulli beam theory is no longer appropriate, and knowledge gained from the Timoshenko beam theory should be used [5, 7-9, 12]. This theory considers concomitantly the bending and shear deflections as well as the mass and the rotatory inertia per unit length.

Timoshenko beam theory will not be developed in detail here. For the purposes of the present work, we will use the theoretical results presented in Ref. [5]. A complete basic text on the theory and some of its applications can be found in Ref. [8]. Additionally, greater elucidation about the theory can be gathered from [7] and [9].

We stress that in the vibration analysis, the concept of slenderness is associated with the distance between null bending moments. Hence, even for slender beams, when the modal order increases, the shear deflection and the rotatory inertia will have an increasing influence on the natural frequencies. We may conclude, therefore, that using the Timoshenko theory is always more suitable when there are short members in the model and/or there is an interest in calculating modes of a higher order. Nevertheless, as shown, the price paid for using this theory reflects the mathematical complexity, which is greater than for the simpler Euler-Bernoulli theory.

Additionally, as the very title of Ref. [5] suggests, it can be interesting to compute the effect of axial load in flexural vibrations. When a beam suffers transverse displacement, its cross-sections suffer translations and rotations in such a manner that the axial loads perform work. If such work is added to the total amount of the elastic energy, the natural frequencies are affected. Therefore, tensile axial loads increase the natural frequencies, while, inversely, compression axial loads decrease the natural frequencies. In this last case, depending on the magnitude of the load, the "dynamic buckling" phenomenon can occur, which is characterized by the total nullification of the bending rigidity of the beam.

The fundamentals of the dynamic problem herein investigated are given in sections 2 and 3, which lead to the DSM formulation given in section 4. Section 5 investigates the numerical approach necessary for the implementation of the secant method [1]. In section 6, this method is modified to accommodate special characteristics of the non-linear transcendental eigenvalue problem explored here.

In this context, two particular numerical aspects are discussed: overflow occurrences and deflation of the DSM determinant.

Because the secant method is an iterative process, its implementation is limited by the computer numerical structure such that overflow occurrences should be predicted and circumvented; without such planning, this process will certainly fail. Therefore, for the calculation of the hyperbolic function arguments as well as for the DSM determinant calculation, the numerical techniques described in section 5.3 should be considered.

For linear eigenvalue problems, whose similarity has been established by the finite elements method [1], the dynamic stiffness matrix determinant is a polynomial of equal or lesser order than the number of degrees of freedom. In this case, polynomial deflation [1, 13] works perfectly as an auxiliary technique of the secant method, allowing for iterative calculation of a certain eigenvalue, in a recursive process, if the lower order eigenvalues are known.

For the present DSM case, nevertheless, its determinant exhibits a non-polynomial form so that the polynomial deflation cannot be successfully utilized. Thus, we introduce in the section 6 the technique called "power deflation": in the calculation of the nth-order eigenvalue, the determinant is deflated by a product of monomial powers; each one of these containing an eigenvalue whose order is lower than n.

2 MEMBER EQUILIBRIUM EQUATIONS

For an elementary volume of length dx (Fig. 1), the equilibrium of lateral forces can be written as:

where F_C (x, t) is the shear force due to bending, m = ρA is the uniform distribution mass and v(x, t) is the transverse displacement.

Analogously, the equilibrium of the elementary volume moments gives:

where M_F (x, t) is the bending moment, j = ρI is the uniformly distributed rotatory inertia, θ(x, t) is the angle of rotation of the cross-section, and P is a static axial load.

The static axial load is supposedly time-invariant, constant throughout the length and positive under tensile condition.

On the other hand, the shearing condition of the elementary volume can be written as:

where G = E/[2(1 + υ)] is the shear modulus and A_S is the shear area that is given as a fraction of the total area of the cross-section (A). The shear area characterizes the stiffness in relation to the deflection due to shear, so that if A_S = ∞ there is no such deflection.

Thus, eliminating F_C (x, t) and M_F (x, t) from these equations, it results in the following differential equation, applicable for both transverse displacement and rotation [5]:

where the dot denotes v(x, t) or θ(x, t). Noted in this equation is the existence of crossed terms that are involved in the relations between the beam parameters: EI = flexural rigidity, GA_S = shear rigidity, m = mass per unit length, j = rotatory inertia and P = static axial load.

Additionally, using the same elementary volume of length dx (Fig. 1), the dynamic equilibrium of axial forces, from which P is excluded because it is constant, can be written as:

where F_N (x, t) is the normal force and u(x, t) is the axial displacement.

3 MEMBER EQUILIBRIUM SOLUTION

In free vibration, the displacements are synchronous and harmonic and can be written as:

stressing that, for the present theory, the rotation θ(x, t) cannot be directly obtained from the derivative of v(x, t), as it could in the Euler-Bernoulli theory.

Substitution of expressions (6) for displacements in the equations from the previous section, after some algebraic transformations, leads to:

where , , , , and are integration constants to be determined by using appropriate boundary conditions, θ_T and θ_H are parameters to be given forward and:

Analogously, in free vibration, the internal forces can be written as:

with

where the parameters F_A ,F_T ,F_H ,M_T and M_H are given forward.

For convenience, we define the following non-dimensional parameters:

where, in Eq. (11) we have:

Thus, we finally write:

where

with the non-dimensional parameters α defined by:

being

4 ASSESSMENT OF THE EXACT DYNAMIC STIFFNESS MATRIX (DSM)

4.1 Element local dynamic stiffness matrix

Employing the convention for the end displacement and the end forces shown in Fig. 1, as well as the origin there indicated for the local reference system, and defining for J = A, T or H :

the dynamic stiffness matrix of an axially loaded Timoshenko member is given by:

where the matrices φ_Q and φ_P are, respectively:

End displacements versus integration constants relationship

End forces versus integration constants relationship

In the expressions (28) and (29), the parameters θ_T , θ_H , F_A , F_T , F_H , M_T and M_H depend on the frequency ω and are defined by Eqs. (18).

By the usual FEM procedures, the matrix of expression (27), after applied to each member of the model, is useful to assemble the classical global DSM [5]. However, for frequencies that coincide with the natural frequencies of the members in the clamped-clamped condition, the classical DSM has poles, not roots.

This can be easily seen by observing that the matrix Φ_Q (ω) has no inverse when ω is one of the natural frequencies of the member in the clamped-clamped condition so that, as consequence, the DSM determinant grows to infinity.

For this reason, the classical DSM is not suitable in an iterative method that needs determinant values to estimate the natural frequencies. To overcome this difficulty, in the following section it is shown how the matrices Φ_P (ω) and Φ_Q (ω) can be used in order to assemble an improved version of the global DSM that has no poles [4]. Then, the improved DSM can be used to extract eigenvalues by the secant method [1], as it will be presented in section 6.

4.2 Improved global dynamic stiffness matrix

By applying the standard FEM procedures, the process by which the matrices φ_P (ω) and φ_Q (ω) are employed for the in-plane model assembly, in the light of the Timoshenko theory, follows the same steps that are described in Ref. [4] for the Euler-Bernoulli theory. For the sake of brevity, the details will not be presented here. The terms are:

M

Total number of members of the model;

N

Total number of nodal degrees of freedom;

Q_o

Vector of the amplitudes of the nodal displacements;

Θ

Vector that contains the constants of integration of all the members;

Matrix that unites the nodal concentrated masses and springs;

ψ _Q

Matrix that determines the connection between the members of the model;

Φ_P

Matrix that unites the matrices φ_P (ω) of all the members;

Φ_Q

Matrix that unites the matrices φ_Q(ω) of all the members.

The eigenvalue problem, which is established through the exact improved dynamic stiffness matrix Ψ(ω), results in the following form [4]:

being necessary to search for ω such that:

For this eigenvalue problem, which does not show poles in the natural frequencies [4], it is possible to obtain the roots of Eq. (31) and then return to Eq. (30) to obtain the corresponding vectors Θ and Q_o. The first of these vectors contains the integration constants, which in the latter instance determine the amplitudes of the internal displacements, and the second vector determines the amplitudes of the nodal displacements.

In these terms, for the calculation of each natural mode shape, the solution of Eq. (30) is based in a method that linearizes the eigenvalue problem in the vicinity of the corresponding natural frequency, as described in Ref. [4]. For the sake of brevity, this method will not be reproduced here.

In Ref. [4], a simple method of determinant tracing was presented to solve Eq. (31). Although robust, it is not computationally efficient because the method requires successive evaluations of the determinant, at frequency intervals not greater than the tolerance specified for the natural frequencies calculations. Thus, the objective is to implement the secant method with the use of power deflation, which is computationally much more efficient. Having in mind this implementation, in the following section some particular issues about the present eigenvalue problem are discussed.

Based on transcendental member stiffness matrices, Ref. [9] defines a new dynamic stiffness matrix, called K_∞ , which also eliminates the poles from the eigenvalue problem. However, this method requires further developments in order to include rigid off-set and end release in its formulation as well as has no direct capability to determine the local mode shapes for which Q_o = 0 and Θ ≠ 0.

5 PARTICULAR ISSUES

As we have shown, derivation of the Timoshenko theory is much more laborious than the Euler-Bernoulli theory; consequently, new physical and mathematical questions arise, which we will discuss in the following sections.

5.1 Limiting compression axial load

By imposing that the denominator in expressions (21) and (22), for the parameters α_T and α_H, respectively, should be real and finite, we may write:

which limits the axial load to:

Thus, the lower the shear rigidity, GA_S , the lower the axial load that can be applied. In practical terms, for actual materials such a limit is sufficiently high to make feasible the calculation for any realistic structure.

On the other hand, even for infinite shear rigidity, the compressive axial load can be physically limited by the buckling phenomenon. As an illustrative example, we may consider the case of a simply supported beam with no shear deflection (A_S = ∞ ⇒ s = 0). For this situation of a Rayleigh beam, the equations in section 3 are applied in the appendix A and give for the natural frequencies:

where, in accordance with the Euler buckling load : P_CR = π²EI/L².

Hence, if the load is compressive (P < 0) and has magnitude greater or equal to the buckling load, null or complex natural frequencies will occur; a fact that characterizes the physical impossibility of dynamic equilibrium.

5.2 Transition frequency

The parameter α_T , present in the arguments of the trigonometric functions, will always result in a real number because the conditions ∂² < 0 and < 0 cannot occur concomitantly.

On the other hand, the parameter α_H , in the arguments of the hyperbolic functions, can result in a real number or a purely imaginary number, without thereby implying any impossibility, either mathematical or physical.

In these terms, a particular frequency value exists for which the transition of α_H occurs from a real number to an imaginary number. Thus, the use of expression (25) and the condition (32) imposes that:

to obtain, from Eq. (11), the following expression for the transition frequency:

The arguments of the hyperbolic functions are purely imaginary numbers for frequencies greater than the transition frequency, a fact that indicates that such functions convert to the trigonometric functions. From a computer calculation point of view, this is a favorable condition because the trigonometric functions, by being limited, do not present overflow problems.

As noted in Eq. (36), if the shear deflection is neglected (s = 0) there is no transition and overflow can occur because the hyperbolic function arguments are always real numbers. The same fact occurs when the rotatory inertia is neglected (r = 0). On the other hand, it can be noted that the transition frequency does not depend on the member length and therefore it is not affected by mesh refinement.

5.3 Treatment of overflow conditions

For computer calculation, there is always an upper limit for the largest number that can be stored in the RAM (R_max ). If during some step of an iterative calculation such a limit is surpassed, this calculation is certainly compromised due to overflow.

The equations in sections 3 and 4 have parameters that increase with the frequency ω such that this may lead to overflow error, particularly because hyperbolic functions are used.

Hence, to overcome overflow occurrences, appendix B is dedicated to the establishment of five possible limiting values (Ω_{lim J} for J = 1,2,3,4,5) for the non-dimensional frequency Ω defined by Eq. (11). Accordingly, each member of the model will have its own limiting frequency given by:

For a framed structure, the overall limiting frequency is the smallest value given by Eq. (37) applied for all the members. If this overall frequency is less than the maximum desired frequency, the only remaining alternative is to promote mesh refinement. As inferred from the example in Fig. 2, the shorter the elements the greater the overall limiting frequency.

 

On the other hand, the more refined the mesh the greater the possibility of overflow occurrence in the calculation of the DSM determinant because the matrix order will be greater. Fortunately, as shown afterwards, scaling operation in the DSM can be employed to avoid such an occurrence.

In order to obtain a more elucidative visual effect, the example in Fig. 2 uses an intentionally low R_max value. In Fig. 2a, it can be observed how each of the non-dimensional limiting frequencies Ω_{lim J} varies with the mesh refinement. By applying Eq. (37), Fig. 2b shows the behavior of the overall limiting frequency ƒ_lim . Finally, Fig. 2c shows how mesh refinement increases the modal order of the greatest natural frequency that can be calculated.

Table 1 explores how the increase in the order of R_max , up to a compatible value with the real computer capacity, affects the modal order of the greatest natural frequency that can be calculated. In general, it is assumed that only unrealistic data should cause relevant computational limitation by overflow error. For example, from the situation in the last line in Table 1, if an unrealistic axial load P = +10¹⁰⁰ ton is applied, the maximal calculable mode falls to an order of 4.

Dynamic stiffness scaling

We consider the dynamic stiffness matrix, which is defined in Eq. (30), now affected by a real and positive scaling constant α_G :

so that the calculation of the DSM determinant can be written as:

where M and N are positive integers, which represent, respectively, the number of elements of the model and the total number of nodal degrees of freedom.

Thus, we note from Eq. (39) that the scaling constant α_G can be used to change the magnitude of the determinant result without, in essence, affecting the eigenvalue problem.

With no changes in the mode shapes and natural frequency results, a scaling constant α_G lower than the unit can avoid the overflow occurrence.

To predict a suitable value of α_G , we apply a process of trial and error, in the following terms:

1. The overflow condition is tested for a value less than R_max , through the appropriate adoption of the safety coefficient S_F > 1: V_max = .

2. In the calculation of the first natural frequency, while the determinant is greater than V_max , the α_G value is diminished and the iterative process is restarted.

3. Once a natural frequency is obtained, the maximum current determinant value (D_max), if greater than V_max , is used to diminish the current scaling constant by using: α_G = α_G [D_max /V_max ]^{1/(6M + N)} .

For some methods other than the secant method, like bisection, scaling of the dynamic stiffness determinant could be conveniently performed during calculation by maintaining its magnitude within an allowable range and storing an additional integer exponent, so avoiding the need of predicting a suitable value for α_G.

6 POWER SECANT METHOD

The objective here is to elaborate a method that serves as an alternative to the tracing method of Ref. [4], leading to a search for the natural frequencies that are computationally more efficient.

We begin with the discussion of the simple secant method; afterwards, we address the question of polynomial deflation, and, finally, we deal with power deflation.

In principle, the simple secant method lends itself to the extraction of only one root of a function. If the function has only one root, the method is directly applied, and it does not need deflation. If, on the other hand, the function has more than one root, the method can still work under the following conditions:

• If the starting points are chosen below the first root, convergence to this root will occur;

• If the starting points are sufficiently close to a given root, convergence to this root will occur.

In the case of functions with several roots, we can apply the deflation technique to extract them in ascending order [1, 13]. Thus, by removing the first root by deflation, if the starting points are chosen from below the second root, convergence to this second root will occur. In succession: being known the ( n-1 ) lower roots, if by deflation these roots are removed, when choosing starting points from below the nth root, convergence to this root will occur.

6.1 Simple Secant Method

Observing the secant method in its simplest form in Fig. 3, we can develop through recursive application of the equation that represents the intersection of the secant line with the λ axis:

Beginning with k = 1, so that it is necessary to arbitrate two starting points λ₁ and λ₂ , the iterations are interrupted when a tolerance test is met, as prescribed by:

where d is the number of significant digits that are desired for the approximation to the sought root. If during the iterations ƒ_k = ƒ_k+1 occurs, the method fails in an unrecoverable way. In general, the method should not be directly applied in case the function has null derivative points, as when more than one root exists. In this context, the deflation technique becomes relevant.

6.2 Polynomial deflation

Lets suppose that the function ƒ(λ), from which we desire to extract the roots, is a polynomial of order N, with N real and distinct roots. Thus, the function is known in an explicit form from its coefficients a_n :

or the function can be imagined in the implicit factorized form, described by monomial products involving its roots

_n:

As shown before, the simple secant method can be applied for calculation of the first root if the starting points that are lower than this root. Once the first root

₁ is obtained, a new function is found by using the following deflation:

Now, the first root of the deflated polynomial F₂ (λ) is in the position of the second root of the original polynomial ƒ(λ). Thus, the secant method can be now applied on the deflated polynomial. If the starting points are chosen below this first root, the method converges to the second root.

Such a deflation procedure can be generalized for the progressive calculation of all the roots.

where F_n (λ) is the deflated form of the original polynomial that makes the secant method to converge to the n-order root, being necessary that all the lower ( n-1) roots have already been calculated.

For the purpose of the eigenvalue extraction of a linear problem, and for the conditions that have been established here, the secant method with polynomial deflation can be applied with success even for the case of repeated natural frequencies [13].

6.3 Power deflation

We have seen in the previous section that the secant method with polynomial deflation works perfectly well for the cases in which the function is a polynomial with real roots. In this context, for linear and non-transcendental eigenvalue problems, similar to the problems that arise from FEM models, the secant method can be employed without reservations. Being K the stiffness matrix and M the mass matrix of a FEM model with N degrees of freedom:

is a polynomial of order N with a maximum of N real and non negative roots [1]. The fact that the polynomial in Eq. (46) is not known in an explicit form, like in Eq. (42), does not hinder the application of the secant method because it is still possible to make numerical evaluations of ƒ(λ). Especially for case of Eq. (46), such numerical evaluations are more efficiently conducted through the Gauss factorization [1]:

where L(λ) is a lower triangular matrix and D(λ) is a diagonal matrix. Once performing factorization, the determinant is simply calculated from the products of the diagonal terms of the matrix D(λ):

On the other hand, for the improved dynamic stiffness matrix Ψ in Eq. (30), the function ƒ(λ) = det{Ψ(λ)} is not a polynomial, although it has non-negative real roots. Thus, the polynomial deflation has little or no possibility of working.

In these terms, alternatively we use the following modified version of Eq. (45), which characterizes what we call power deflation:

where each N_j power should be determined in a manner that, given two starting points λ_1n and λ_2n , the following condition is met:

For calculation of the first eigenvalue (n = 1 ⇒ j = 0), we have F₀ (λ) = f(λ) and _o = 0. Consequently, before starting the iterations for each root, it is guaranteed that the first secant line favorably intersects the λ axis at a point located ahead the starting points. Then, from Eq. (50), we can obtain the following general expression for calculating the N_j powers that are applied in Eq. (49):

If this expression results in N_j < 1, N_j = 1 is adopted because the corresponding deflation is not needed.

For the iterations of the first mode (n = 1⇒ j = 0), we calculate N₀ using Eq. (51) and apply Eq. (49) in order to obtain F₁ to use in the secant method. For the iterations of the second mode (n = 2 ⇒ j = 1), using F₁, we calculate N₁ by Eq. (51) and apply Eq. (49) in order to obtain F₂ to use in the secant method; and thus, the other following modes are progressively calculated.

To illustrate the operation of the secant method with power deflation, or here simply called "the power secant method", we present in Fig. 4 the results for a simple clamped-clamped Euler-Bernoulli beam. In this case, we used ζ = 2 and appropriate starting points in order that convergence occurs for up to ten iterations. In the figure, we can visualize how the deflated functions converge for the eigenvalues of the problem, which are localized at the peaks of the function -log{|ƒ(λ)|}.

Obviously, the success of the method, in terms of the number of iterations that are needed to achieve convergence, strongly depends on the starting points and on the appropriate choice of the parameter ζ.

Depending on the problem, the parameter ζ can be previously admitted as a constant or varied, by trial and error, until the iterative process demonstrates appropriate convergence. In this last case, in the beginning of the iterations for each root, ζ should be re-estimated whenever necessary.

In Table 2, the method is applied to the ten possible distinct cases of boundary conditions for a simple Euler-Bernoulli beam, for which a total of six different frequency equations exist. Without the benefit of prior knowledge of the roots, the solution for each one of these equations uses arbitrary starting points as well as automatically determined value for the parameter ζ.

We note in Table 2 that, under this more general situation useful for computer programming, there is a reasonable increase in the number of iterations. We also observe that the final results precisely coincide with the known solutions of the problems [2].

In spite of the number of iterations, the power secant method is always faster than the bisection method. For comparison, lets take the example of the clamped-clamped beam in Table 2. The power secant method performed 168 evaluations of the frequency equation while a bisection method needs more than 1000 evaluations to reach the same accuracy.

7 COMPARISON EXAMPLES

The power secant method was implemented in the VIGENE program [3], which was developed on the MATLAB platform. This program is a part of the VIANDI computational system that is available for download at http://www.gmsie.usp.br.

Section 7.1 presents a very simple example, which has analytical solution. Subsequently, sections 7.2 and 7.3 explore other simple examples, with results from the literature [7, 8, 12] also compared with additional solutions by FEM where, by dividing the members into hundred equal parts, the mesh was deemed appropriate for the comparisons.

Finally, in order to explore the entire VIGENE capability, section 7.4 presents a more complex model.

7.1 Simply supported Euler-Bernoulli beam under axial load

As shown in the appendix A, if shear deflection and rotatory inertia are disregarded, the simply supported beam has exact analytical solution, which does not require application of any iterative numerical processes. We can observe in Table 3 the perfect concordance between this analytical solution and the numerical solution of the VIGENE program.

7.2 Single Timoshenko beam for various boundary conditions

For comparison, from [7] and [12] we took the examples shown in Tables 4 and 5, respectively. A good concordance can be observed between VIGENE and ANSYS programs.

7.3 Building like frames comprising Timoshenko beam members

For the examples in Fig. 5 and 6 from [8], all the members have the same material and section properties, as indicated in Tables 6 and 7, respectively. An excellent agreement for the results of both examples can be observed, even when compared with FEM solutions by ANSYS.

7.4 Robotic handler comprising Timoshenko beam members

In order to explore all the potential of the present work, Fig. 7 depicts an actual example of a robotic handler where effects like rigid offset, end force release and skewed edge are employed. Table 8 gives section and material properties for all the members while Table 9 shows the obtained natural frequencies. The results of both VIGENE and ANSYS programs agree very well but it must be pointed out that the ANSYS model needed more than 240 beam elements to achieve the same accuracy obtained by the exact solution of the program VIGENE, which uses only four elements.

8 CONCLUSIONS

This work consolidates the possibility of using an iterative process for extracting eigenvalues, without the compulsory use of the Wittrick-Williams algorithm [6, 10, 11, 14]. Ref.[4] gives more details of the reasons why such an algorithm has operational difficulties when dealing with models that incorporate the special effects of end release and rigid offset.

For the examples here presented, the results in Table 10 show how an iterative determinant search method, based on the improved DSM of section 4.2, can be more computationally efficient. In terms of CPU time, the power secant method is always faster than a bracketing method based on the Wittrick-Williams algorithm. Therefore, as a well-known fact [1], the secant method is faster than any suitable determinant tracing, bracketing or bisection method.

For CPU time comparison with FEM solutions, it must be considered that the FEM models, in general, cannot give accurate results unless fine meshes are employed. In the example of Fig 7, for instance, more than 240 beam elements were necessary in order to achieve a good accuracy. Having this in mind, Table 10 shows that the power secant method can be faster than the subspace iteration method [1] employed by the ANSYS program. When compared with FEM, however, the main advantage of the power secant method based on the improved DSM is the capability to find exact solutions that do not depend on the mesh refinement.

Finally, it must be emphasized that the Wittrick-Williams algorithm must use the classical formulation for the DSM while the power secant method must use the improved formulation as defined in section 4.2. As an additional advantage, the improved DSM preserves any local mode shape for which the nodal displacement vector Q_o is null. Any other method based on the classical DSM cannot appropriately detect local modes since the classical DSM becomes ill-conditioned. Moreover, using an inverse vector iteration method based on the improved DSM, it is possible to detect all the distinct mode shapes even when repeated natural frequencies exist.

Received 23 Mar 2010;

In revised form 21 Jun 2010

APPENDIX A. SIMPLY SUPPORTED RAYLEIGH BEAM UNDER AXIAL LOAD

In the case of disregarding shear deflection (s = 0 ⇒ ξ = 1 ; η = r and f_T = ∞), the problem possesses an analytical solution that does not require the application of iterative numerical processes. By applying the boundary conditions in the functions that regulate the transverse displacements (7) and bending moments (10), respectively given by:

we have :

a) Null transverse displacement at x = 0

b) Null bending moment at x = 0

such that from Eq.s (A.3) and (A.4) it results that: = = 0.

c) Null transverse displacement at x = L

d) Null bending moment at x = L

As the transition frequency is infinity, we have S_H (L) ≠ 0. Hence, Eq.s (A.5) and (A.6) are simultaneously satisfied if = 0 and S_T (L) = 0, being ≠ 0 and arbitrary. Thus, the non-trivial solutions are obtained from:

It can be observed that the result of Eq. (A.7) is independent of the parameters r, s and p, respectively defined by Eqs. (12), (13) and (14), in a manner that it is also valid for any case of simply supported Timoshenko beam. The difference for the simply supported Euler-Bernoulli beam, with no axial loading, resides in the way by which the eigenvalues α_T = iπ lead to the natural frequencies. For an unloaded Euler-Bernoulli beam, with r = s = p = 0, the solution simplifies to Ω = (iπ)².

In order to use the known α_T , however, the Rayleigh beam (p ≠ 0 and r ≠ 0) requires some elaboration of the solution for the natural frequency. Then, from Eqs. (24) and (25), with ξ = 1 and η = r, it follows that δ² = 4Ω² and = p² +r²Ω² . So that, using (11) and (13), after the following deduction:

for the non-dimensional frequency we may write :

Thus, using α_T = iπ in Eq. (A.8), as well as the expressions (11), (14) and (17), the natural frequencies are given by:

By converting Eq. (A.9) for frequencies in Hz and by using the Euler buckling load (P_CR = π²EI/L² ), we can finally write:

APPENDIX B. OVERFLOW LIMITING FREQUENCIES

B.1 Limiting frequency (Ω_{lim 1}) due to

Observing Eq. (24), the condition:

results in the following limit for the non-dimensional frequency:

B.2 Limiting frequencies (Ω_{lim 2} , Ω_{lim 3}and Ω_{lim 4} ) due to

The expression (23) for can be viewed in the following form:

where

From this, three possible limits for the non-dimensional frequency can occur: 1.

B.3 Limiting frequency (Ω_{lim 5}) due to hyperbolic functions

If the current frequency is greater or equal to transition one (ƒ > ƒ_T ), overflow in the hyperbolic functions will not occur because α_H , according to Eq. (22), is purely imaginary number. In contrast, (ƒ < ƒ_T ⇒ δ² > 0 ⇒ > ) α_H is real and must be lower than the following computational limit:

Thus, using Eqs. (22), (24), (B.4) and (B.5), from the condition α_H< α_Lim , we arrive at:

from where the following possible cases give a new limiting value for the non-dimensional frequency (Ω(ω) < Ω_{lim 5}):

Case 1 If r² > 0 and s² = 0

Case 2 If r² = 0 and s² > 0

Case 3 If r² > 0 and s² > 0 so

Case 4 If r² = 0 and s² = 0 so

For expressions (B.11), (B.13) and (B.15), overflow does not occur because the hyperbolic functions are such that the limiting value of the non-dimensional frequency results in a number that is equal to the maximum number that can be stored.

• [1] K. Bathe. Finite Element Procedures Printice-Hall, New Jersey, 1996.
• [2] R.D. Blevins. Formulas for Mode Shape & Natural Frequency Krieger, Malabar, Florida, 1979.
• [3] C.A.N. Dias. Computational system for dynamic analysis of in-plane beam structures published on line in http://www.gmsie.usp.br
• [4] C.A.N. Dias and M. Alves. A contribution to the exact modal solution of in-plane beam structures. Journal of Sound and Vibration, 328(4-5):586606, 2009.
• [5] W.P. Howson and F.W. Williams. Natural frequencies of frames with axially loaded timoshenko members. Journal of Sound and Vibration, 26(4):503515, 1973.
• [6] S. Ilanko and F.W. Williams. Wittrick-Williams algorithm proof of bracketing and convergence theorems for eigenvalues of constrained structures with positive and negative penalty parameters. International Journal for Numerical Methods in Engineering, 75:83102, 2008.
• [7] L. Majkut. Free and forced vibrations of timoshenko beams described by single difference equation. Journal of Theoretical and Applies Mechanics, 47(1):193210, 2009.
• [8] J.F. Martins. Rotatory inertia and shear force effects on the natural frequencies and dynamic response of framed structures PhD thesis, Engineering School of São Carlos, University of São Paulo, Brazil, 1998.
• [9] F.W. Williams, D. Kennedy, and M.S. Djoudi. Exact determinant for infinite order FEM representation of a Timoshenko beam-column via improved transcendental member stiffness matrices. International Journal of Numerical Method in Engineering, 59:13551371, 2004.
• [10] F.W. Williams, S. Yuan, K. Ye, D. Kennedy, and M.S. Djoudi. Towards deep and simple understanding of the transcendental eigenproblem of structural vibrations. Journal of Sound and Vibration, 256(4):681693, 2002.
• [11] W.H. Wittrick and F.W. Williams. A general algorithm for computing natural frequencies of elastic structures. Quarterly Journal of Mechanics and Applied Mathematics, 24:263284, 1971.
• [12] J-S. Wu and D-W. Chen. Free vibration analysis of a Timoshenko beam carrying multiple spring-mass systems by using numerical assembly technique. International Journal of Numerical Method in Engineering, 50:10391058, 2001.
• [13] C-D. Yuan and W-H. Chieng. Method for finding multiple roots of polynomials. Computers and Mathematics with Applications, 51:605620, 2006.
• [14] S. Yuan, K. Ye, and F.W. Williams. Second order mode-finding method in dynamic stiffness matrix method. Journal of Sound and Vibration, 269:589708, 2004.

Publication Dates

• Publication in this collection
  06 Mar 2012
• Date of issue
  Sept 2010

History

• Received
  23 Mar 2010
• Accepted
  21 June 2010