Services on Demand
Journal
Article
Indicators
Related links
 Cited by Google
 Similars in SciELO
 Similars in Google
Share
Latin American Journal of Solids and Structures
Online version ISSN 16797825
Lat. Am. j. solids struct. vol.7 no.3 Rio de Janeiro Sept. 2010
http://dx.doi.org/10.1590/S167978252010000300005
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 – 05508900 – 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 illconditioned; 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 WittrickWilliams 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 EulerBernoulli skeleton structures is discussed at length by Dias and Alves [4]. According to the EulerBernoulli 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 EulerBernoulli beam theory is no longer appropriate, and knowledge gained from the Timoshenko beam theory should be used [5, 79, 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 EulerBernoulli 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 crosssections 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 nonlinear 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 nonpolynomial 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 nthorder 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 crosssection, and P is a static axial load.
The static axial load is supposedly timeinvariant, 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 crosssection (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 EulerBernoulli 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 nondimensional parameters:
where, in Eq. (11) we have:
Thus, we finally write:
where
with the nondimensional 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 clampedclamped 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 clampedclamped 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 inplane model assembly, in the light of the Timoshenko theory, follows the same steps that are described in Ref. [4] for the EulerBernoulli 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 offset 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 EulerBernoulli 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} = π^{2} EI/L^{2}.
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 ∂^{2} < 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 nondimensional 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 nondimensional 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^{100} 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 ( n1 ) 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 λ_{1} and λ_{2} , 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 _{1} is obtained, a new function is found by using the following deflation:
Now, the first root of the deflated polynomial F_{2} (λ) 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 norder root, being necessary that all the lower ( n1) 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 nontranscendental 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 nonnegative 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_{0} (λ) = 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_{0} using Eq. (51) and apply Eq. (49) in order to obtain F_{1} to use in the secant method. For the iterations of the second mode (n = 2 ⇒ j = 1), using F_{1}, we calculate N_{1} by Eq. (51) and apply Eq. (49) in order to obtain F_{2} 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 clampedclamped EulerBernoulli 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 reestimated whenever necessary.
In Table 2, the method is applied to the ten possible distinct cases of boundary conditions for a simple EulerBernoulli 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 clampedclamped 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 EulerBernoulli 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 WittrickWilliams 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 WittrickWilliams algorithm. Therefore, as a wellknown 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 WittrickWilliams 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 illconditioned. 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.
References
[1] K. Bathe. Finite Element Procedures. PrinticeHall, New Jersey, 1996. [ Links ]
[2] R.D. Blevins. Formulas for Mode Shape & Natural Frequency. Krieger, Malabar, Florida, 1979. [ Links ]
[3] C.A.N. Dias. Computational system for dynamic analysis of inplane beam structures. published on line in http://www.gmsie.usp.br. [ Links ]
[4] C.A.N. Dias and M. Alves. A contribution to the exact modal solution of inplane beam structures. Journal of Sound and Vibration, 328(45):586–606, 2009. [ Links ]
[5] W.P. Howson and F.W. Williams. Natural frequencies of frames with axially loaded timoshenko members. Journal of Sound and Vibration, 26(4):503–515, 1973. [ Links ]
[6] S. Ilanko and F.W. Williams. WittrickWilliams 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:83–102, 2008. [ Links ]
[7] L. Majkut. Free and forced vibrations of timoshenko beams described by single difference equation. Journal of Theoretical and Applies Mechanics, 47(1):193–210, 2009. [ Links ]
[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. [ Links ]
[9] F.W. Williams, D. Kennedy, and M.S. Djoudi. Exact determinant for infinite order FEM representation of a Timoshenko beamcolumn via improved transcendental member stiffness matrices. International Journal of Numerical Method in Engineering, 59:1355–1371, 2004. [ Links ]
[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):681–693, 2002. [ Links ]
[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:263–284, 1971. [ Links ]
[12] JS. Wu and DW. Chen. Free vibration analysis of a Timoshenko beam carrying multiple springmass systems by using numerical assembly technique. International Journal of Numerical Method in Engineering, 50:1039–1058, 2001. [ Links ]
[13] CD. Yuan and WH. Chieng. Method for finding multiple roots of polynomials. Computers and Mathematics with Applications, 51:605–620, 2006. [ Links ]
[14] S. Yuan, K. Ye, and F.W. Williams. Second order modefinding method in dynamic stiffness matrix method. Journal of Sound and Vibration, 269:589–708, 2004. [ Links ]
Received 23 Mar 2010;
In revised form 21 Jun 2010
* Author email: candias@usp.br
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 nontrivial 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 EulerBernoulli beam, with no axial loading, resides in the way by which the eigenvalues α_{T} = iπ lead to the natural frequencies. For an unloaded EulerBernoulli beam, with r = s = p = 0, the solution simplifies to Ω = (iπ)^{2}.
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 δ^{2} = 4Ω^{2} and = p^{2} +r^{2} Ω^{2} . So that, using (11) and (13), after the following deduction:
for the nondimensional 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} = π^{2} EI/L^{2} ), 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 nondimensional 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 nondimensional 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} ⇒ δ^{2} > 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 nondimensional frequency (Ω(ω) < Ω_{lim 5}):
Case 1 If r^{2} > 0 and s^{2} = 0
Case 2 If r^{2} = 0 and s^{2} > 0
Case 3 If r^{2} > 0 and s^{2} > 0 so
Case 4 If r^{2} = 0 and s^{2} = 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 nondimensional frequency results in a number that is equal to the maximum number that can be stored.