General exact harmonic analysis of In-Plane Timoshenko Beam Structures

The exact solution for the problem of damped, steady state response, of in-plane Timoshenko frames subjected to harmonically time varying external forces is here described. The solution is obtained by using the classical dynamic stiffness matrix (DSM), which is non-linear and transcendental in respect to the excitation frequency, and by performing the harmonic analysis using the Laplace transform. As an original contribution, the partial differential coupled governing equations, combining displacements and forces, are directly subjected to Laplace transforms, leading to the member DSM and to the equivalent load vector formulations. Additionally, the members may have rigid bodies attached at any of their ends where, optionally, internal forces can be released. The member matrices are then used to establish the global matrices that represent the dynamic equilibrium of the overall framed structure, preserving close similarity to the finite element method. Several application examples prove the certainty of the proposed method by comparing the model results with the ones available in the literature or with finite element analyses.


) (x F 1 INTRODUCTION
Many modern structures are formed by beam elements.These skeleton like structures are subjected to static and dynamic loads.Their beam members can be of various sizes, including beams with small length to beam height ratio.For the analysis of these structures, it is important to use a more refined beam theory, where the assumption of the cross section to remain plane is not enforced.Besides, harmonic loads can be of high frequency, when then it is important to keep in the beam model the cooperation of the rotatory inertia to the overall structure response.This motivates the use of the Timoshenko beam theory to obtain the damped steady state response for general plane frames subjected to harmonically time varying external forces. The study reported here concerns with an exact harmonic analysis using the classical Dynamic Stiffness Matrix, DSM.The problem at hand is non-linear and transcendental with respect to the excitation frequency [Howson and Williams (1973)].The approach used to solve it is by the use of the Laplace transform.
Focusing on the calculation of natural frequencies, Howson and Williams (2003) present a formulation based on the classical DSM obtained by a set of decoupled fourth order partial differential equations (PDE) for the unknowns deflection and rotation of the beam cross section.Dias and Alves (2009) also derive the DSM via the same procedure but reaches an improved formulation, argued to be more suitable for the eigenproblem solution.It has been noticed [Schanz and Antes (2002)] that the dynamic analyses of beams can be performed by decoupling deflection and rotation.In the study described here, the original coupled PDEs, combining deflection, rotation, bending moment and shear force, are directly employed in order to reach the DSM formulation as well as the formulation for the equivalent load vector that arises due to distributed loads on the beams and frames.
We remark that, to be considered exact, a solution must adopt no assumptions other than those of the beam theory itself.Hence, if the mode superposition method is employed, the result is not exact once it is affected by series truncation error.For this reason, many good publications Latin American Journal of Solids and Structures 11 (2014) 2171-2202 dedicated to dynamic frame response calculation do not conform to the requirements of an exact solution and therefore are not considered here.Concerning the dynamic response of frames, Table 1 compares different solution methods found in the literature [Abu-Hilal (2003), Foda and Albassam (2006), Gürgöze and Erol (2001), Lin (2008), Loudini et. al. (2006), Majukt (2009), Mehri et. al. (2009), Saeid (2011), Schanz and Antes (2002), Tang (2008)].In the present context, it is imperative that the PDE representing the governing equations of a given in-plane structure subjected to harmonic forces has to be solved exactly.This can be achieved by using Laplace transform [Abu-Hilal (2003), Foda andAlbassam (2006), Loudini et. al. (2006), Saeid (2011)] and/or Green functions [Abu-Hilal (2003), Tang (2008), Davar and Rahmani (2009)].Gürgöze and Erol (2001) use a distinct method based on a receptance matrix but cannot be considered exact since series solution is employed.
As seen in Table 1, only the solutions given in (Shanz and Antes, 2002;Foda and Albassam, 2006;Lin, 2008;Majkut, 2009;Saeid, 2001) could be generalized in order to solve models of arbitrary numbers of beams and boundary conditions.If further features are to be considered, e.g.concentrated masses and springs, these solutions would be limited to the ones developed in (Foda and Albassam, 2006;Majkut, 2009).Even so, these references are dedicated to solve single and/or in-line beams.When considering the capability to solve framed structures with the features of concentrated masses and springs and rigid bodies attached to the members ends, only Seeid (2001) can be highlighted.These features are not taken into account by Antes et al. (2004), which also deals with harmonic loads applied to Timoshenko frames.Mei (2008Mei ( , 2012) ) present an interesting model that considers in-plane vibration of some restricted forms of frames.The analyses presented in these references are developed, as here, along the Timoshenko bending theory.In contrast, the present paper is more general inasmuch as it handles any shape of portal planar frame subjected to harmonic loads and it solves the equations of motion via the Laplace transform.
Except for the case of transient analysis, which is beyond the scope of this study, the solution of the present investigation is unique in the sense that attends to all of the requisites listed in Table 1.Indeed, to best of the author knowledge, no other publication in the context of harmonic analysis pays attention to effects like rigid offset and end release.Few publications (Abu-Hilal, 2003;Tang, 2008;Saeid, 2001) have considered distributed loads and even fewer have incorporated damping (Loudini, 2006;Abu-Hilal, 2003) in their analysis.
The framed structure under consideration is composed by in-plane members, which are connected to each other through previously defined nodes.Each member has a flexible span, which obeys the Timoshenko beam theory, and may also have rigid bodies attached to its ends, where forces can be released.Fig. 1 depicts the internal forces and displacements for the flexible span, while Fig. 2 shows the end forces, end displacements, distributed loads, and attached rigid bodies for a typical member.Since the rigid bodies may or may not have mass, herein they are alternatively named rigid offsets.Taking into account external and internal damping effects, Section 2 is dedicated to presenting the basic equilibrium equations of a flexible span, in both axial and transverse directions.For the transverse equilibrium, the Timoshenko beam theory (Howson and Williams, 1973;Dias and Alves, 2009;Schanz and Antes, 2002;Loudini et al., 2006) is applied in order to consider the effects of shear deflection and distributed rotatory inertia.In Section 2.2, under the assumption in Section 2.1 that all the excitations vary harmonically with the same frequency, the steady state response is formulated in terms of displacement and force amplitudes, which are then subjected to the Laplace transform.Although this assumption reduces the application scope, it characterizes the well-known harmonic analysis, useful to identify harmful vibrations due to resonance occurrences.On the other hand, considering that a given beam like structure is linear, the proposed method can be extended to a more general excitation type, with more than one distinct excitation frequency.The first order governing equations of Section 2 are a coupled PDE system that might be manipulated in order to obtain uncoupled PDEs for the displacements (Howson and Williams, 1973;Foda and Albassam, 2006).However, this operation would produce fourth order PDEs that are less suitable for Laplace transform (Schanz and Antes, 2002).Hence, in Section 3, the equations of Section 2 are directly subjected to Laplace transforms such that, after some easier but laborious algebraic operations, decoupled results for displacements and forces could be obtained in the Laplace domain.
After applying inverse Laplace transforms to the results of Section 3, Section 4 gives the desirable formulations for the internal displacements and forces of the member flexible span, which are expanded to the model global matrices in Section 5.As it is shown there, from these formulations it is possible to obtain the relation between end forces and end displacements such that the member dynamic stiffness matrix and the member equivalent load vector can be defined.Then, using a similar technique to the finite element method, it is shown how these member matrices are employed to establish the global matrices that represent the dynamic equilibrium of the framed structure.
Latin American Journal of Solids and Structures 11 (2014) 2171-2202 The solution of this equilibrium system directly gives the nodal displacements, from which easily follows member stresses.Some careful chosen examples in Section 6 illustrate the major features of the present method.

BASIC EQUATIONS
According to the Timoshenko beam theory (Howson and Williams, 1973;Dias and Alves, 2009;Schanz and Antes, 2002;Loudini et al., 2006) and considering the internal forces and displacements in Fig. 1, with the distributed linear varying external forces in Fig. 2.b, the following equations are obtained, with the prime denoting the derivative wrt position, x , overdot denoting derivative wrt time, t , and all the loads varying in time with the same circular frequency  : Latin American Journal of Solids and Structures 11 (2014) 2171-2202 Normal stress, σ N ( , ) [ ( , ) ( , )] Transverse force equilibrium Here, E and G are the elastic and shear moduli, u and v are the axial and total transverse beam displacement, v s is the transverse displacement due to shear, A and A s are the total cross-section and shear areas, a A and a L are the slope of the distributed axial and transverse load, b A and b L are the distributed axial and transverse load, m is the mass per unit length, p and P are distributed axial loads.
Latin American Journal of Solids and Structures 11 (2014) 2171-2202 In these equations, the internal and external damping coefficients, I c and E c , respectively, are defined in order to reproduce the Rayleigh formulation, where E c represents the mass proportional coefficient and I c the stiffness proportional coefficient.From this, by using the orthogonal condition of the natural modes, it is possible to establish the following expression for the modal damping ratio i Meaningful values for I c and E c can be obtained whenever a pair of modal damping ratio values is known.Therefore, after applying Eq. ( 10), these damping coefficients can be computed by: where the natural frequencies ( On the other hand, by imposing that both Rayleigh damping coefficients must be non-negative and that the peak frequencies pi  given by: must be real, the following restrictions must be obeyed: The external damping coefficient E c accounts for the environment viscous actions expressed by the forces The internal damping coefficient I c accounts for the energy dissipation due to the structure deformation rate and can be related with the Kelvin-Voight damping coefficient V K (Shanz and Antes, 2002): A more complete characterization of damping in dynamic structural systems can be found in Oliveto, and Greco (2002), where it is shown that the Rayleigh coefficients are independent of the boundary conditions, no matter whether the system is continuous or discrete.

Steady state response
For the steady state response, all the previous defined internal forces and displacements can be written according to the general equation with time and space being decoupled.Performing time derivatives in equations (1-9), the amplitudes are related by: where, for instance, ( ) v x , is the amplitude of the beam total deflection.

Laplace Transform
Applying the Laplace transform over the previous defined amplitudes: the differential equations ( 16) are substituted by the following simple algebraic expressions: It is adopted in the sequence, according to Fig. 2.a, the convention for the end forces and end displacements as Latin American Journal of Solids and Structures 11 (2014) 2171-2202 Now, solving Eqs.(18.a) and (18.b) for the axial displacement results in: Total deflection can be obtained after considerable algebraic manipulation of Eqs.(18.c-f), yielding: Latin American Journal of Solids and Structures 11 (2014) 2171-2202 Analogously, solving Eqs.16(c-f) for the shear force and bending rotation results in: The Laplace transform of the bending moment can be obtained by inserting Eq. ( 28 A common practice in the specialized literature is to decouple Eqs. ( 16) in order to obtain a set of differential equations involving only one unknown function.For the Timoshenko theory, complicated terms involving fourth order derivatives will appear, so the use of the Laplace transform is troublesome.To avoid this, here the Laplace transform was directly applied to the coupled first order Eqs.( 16) by applying the transform to second order coupled equations involving only deflection and rotation.As a remark, Schanz and Antes (2002) also avoids to work with fourth order derivatives.

INVERSE LAPLACE TRANSFORM
By applying the inverse transform to Eq. ( 20) for the axial displacement and then using the result in Eq. (16.a) for the normal force, these corresponding amplitudes are then given, respectively, by: where, knowing that As for the total deflection, it can be obtained by applying the inverse transform to Eq. ( 23), yielding The rotation due to bending is obtained by substituting Eq. ( 23) into Eq.( 28) and applying the inverse transform so that the amplitude is: Latin American Journal of Solids and Structures 11 (2014) 2171-2202 Finally, the shear force amplitudes come from the inverse transform of Eq. ( 27), resulting in while the bending moment amplitude comes from the use of Eq. (37) in expression (16.e),

Flexible span equilibrium
Applying the following boundary conditions for the end displacement (see Fig. 2.a): after laborious algebraic work, the equations of the previous section can be rearranged in the following expressions for the internal displacements: with the functions on the right hand sides (from Analogously, applying the following boundary conditions for the end forces: the internal forces are given by: with the functions on the right hand sides (from and using the previous equations (43.a) for the axial displacement and (45.a) for the normal force, the axial equilibrium is expressed by: 48), after multiplication by the inverse of PA φ , gives the following classical matrix equilibrium equation:

Now, defining the vectors:
with the vector of equivalent axial loads being while the axial dynamic stiffness matrix is Analogously, transverse equilibrium can be written as with the equivalent lateral load vector and the dynamic stiffness matrix given by, respectively,

With
Latin American Journal of Solids and Structures 11 (2014) 2171-2202 Now, combining both axial and lateral matrix equilibrium equations, it follows that: with the end displacements and end forces vectors being given, respectively, by The equivalent load vector is and the dynamic stiffness matrix is Latin American Journal of Solids and Structures 11 (2014) 2171-2202

Rigid offset
Observing Fig. 2, the conversions from the flexible span ends to the member ends are given by: Now, collecting P ˆ and Q ˆ from expressions (64, 65) and substituting into Eq.( 59), member equilibrium requires that: where, for the member equivalent load vector and the member dynamic stiffness matrix, respectively, we have: Latin American Journal of Solids and Structures 11 (2014) 2171-2202

End force release
Let P I denote a permutation matrix that moves, to the last line, the degree of freedom to be released.If, for any end displacement vector Q , it is imposed that the corresponding end force must be null, it can be written that: By doing so, the released versions for the member equivalent load vector and for the member dynamic stiffness matrix are, respectively: with the partitions Ea P ~, Finally, the displacement of the released degree of freedom can be obtained from:

MODEL GLOBAL MATRICES
The member equilibrium equations seen in the previous sections can be also expressed as with the subscript i added to identify the i-th member of the model and the dependency on the excitation frequency  is explicitly indicated.
Then, the usual FEM model assembly leads to the global dynamic equilibrium written as: where ) ( P is the vector of the applied nodal loads, ) ( Q is the vector of the nodal displacements, both in global coordinates, and is the global equivalent load vector, with i R denoting the member rotation matrix and i ψ denoting the connection matrix that relates the member degrees of freedom with the model degrees of freedom.The global dynamic stiffness matrix is expressed by: with the backslash operator denoting an appropriate solution method, based, for instance, on the Gauss elimination.Once the vector ) ( Q has been obtained, the member end displacements are evaluated by: and the member end forces by: Latin American Journal of Solids and Structures 11 (2014) 2171-2202

APPLICATION EXAMPLES
The exact harmonic analysis proposed herein was implemented in the standalone application VIHAND, developed in the MATLAB platform.This code is part of the VIANDI computational system for the exact static and dynamic analysis of skeleton-like structures (Dias, 2011) and freely available at www.gmsie.usp.br.In order to explore the above exact solution, some examples were judiciously chosen and compared with other analytical or numerical approaches from the literature, as now described.

Euler-Bernoulli beam
When shear deflection and rotatory inertia are disregarded, the Timoshenko theory must reproduce the Euler-Bernoulli theory.Therefore, the aim of this section is to verify the present method for the particular case when The examples given in Tables 2 and 3, where no damping is presented, show a very good agreement with the exact solution given in (Abu-Hilal, 2003), which uses Green functions obtained by the Laplace transform.In both tables, additional results obtained by the approximate mode superposition method for different number of modes are also listed.Clearly, it can be seen that the mode superposition series solution improves with increasing number of modes and that the convergence is quite good for the deflection but not so good for the shear force.This confirms the well-known fact that the mode superposition method lacks accuracy for predicting variables involving derivatives wrt position.For the damped cases in Table 4, Abu-Hilal (2003) does not supply written tabular results.For this reason, the results were obtained by using the formulae presented there.Additionally, a FEM model was employed to solve the problem.1000 beam elements have to be used to model the beam so to reach the same accuracy as that given by the present method.As can be observed in the table, the agreement of the three methods is quite good.
Based on a receptance matrix, Gürgöze and Erol (2001) offer a series solution for the example in Table 5 as well as an exact solution based on boundary value formulation.Again, the agreement with the present method is very good.

Timoshenko beam
For the examples with no damping listed in Table 6, Foda and Albassam (2006) use Green functions to provide exact solutions for Timoshenko beams with attached masses and springs.Unfortunately, the authors do not supply written tabular results, so only approximations could be collected from the graphics presented there.As seen in the  By adding damping effects, the same problem has the frequency response depicted in Fig. 3, where a very good agreement with the FEM solution can be observed.

Elastic robot manipulator
The authors did not find in the open literature examples of the exact behavior of more complex skeleton structures.They would allow exploring further the resources in VIHAND and the theory here presented.However, a more refined model representing a robot manipulator is suggested in Fig. 4, whose free vibration response is studied in (Dias, 2010).Here, the effects of rigid offset, bending release, skewed edge, static and dynamic loads, concentrated mass, long and short beams and an external spring, besides damping, are all present.The results obtained by the theory in this article are then compared with the FEM, which is the only alternative that permits to solve such a problem.It is interesting to notice that VIHAND is capable of dealing with the interaction between static and dynamic loads.In the case of this robot structure, when the horizontal arm of the manipulator receives a static compression load, the natural frequencies of the structure decreases.This leads to a harmful vibration due to the resonance between the actuator excitation frequency and the first natural mode.As can be seen in Fig. 5, between 12 and 13 Hz the dynamic stress increases more than twofold, a fact that can severely compromise the structure fatigue life.The results of the VIHAND and ANSYS programs agree very well.However, it must be pointed out that the ANSYS model needed more than 240 beam elements to achieve the same accuracy as that obtained by the exact solution of the program VIHAND, which employs only four elements.Two aspects of the present model are worth considering.First, the method is not suitable for spatial frames analyses, at least in an exact way.The reason is that, for these structures, it is likely that torsion-flexure coupling will be important and torsion is not here taken into account.The consideration of this coupling is far from simple and no attempt was made in the present study to handle it.The second point to highlight is that this model does not intend to be an alternative to the finite element method.But, being possible to obtain an exact solution, it is certain that finite element procedures can benefit from this study inasmuch as the solution procedure here outlined can serve as a reference for this and other numerical procedures.

CONCLUSIONS
By applying a Laplace transform to the coupled first order PDEs, which combine displacement and force amplitudes, this work consolidates an exact solution method for the frequency response of in-plane framed structures.The method, based on the Timoshenko beam theory, is quite general since it allows the use of concentrated masses and springs, rigid bodies, end force release, skewed edges, internal damping and external damping as well as concentrated and distributed harmonic loads.
The reliability of the proposed method could be tested by comparing it against simple beam examples found in the literature.The method was also fully verified by means of a more complex model involving the majority of its capabilities and comparing the results with the FEM.
Although a priori developed for a single excitation frequency, the method can be easily extended for more complex and arbitrary time varying excitations by using Fourier series and/or Fourier transform.Hence, the method is useful for spectral analysis implementation in the context of random vibrations.

Figure 1 :
Figure 1: Timoshenko beam: internal forces and displacements at local coordinates.

K
with C denoting the matrix of the nodal concentrated springs and C M denoting the matrix of the nodal concentrated masses.Once the global equilibrium equation (80) has been assembled, the solution is obtained by:

Figure 3 :
Figure 3: Frequency response for the simply supported Timoshenko beam of Table 6 with two symmetrically attached springs: (a) deflection at load location [mm]; (b) bending moment at load location [Nm].Damping: 10 .0 1  

Figure 5 :
Figure 5: Frequency response for the elastic robot manipulator of Fig. 4: (a) deflection [mm] at point B; (b) bending moment [ton.cm] at point A and (c) shear force [ton] at point A.

Table 1 :
Comparison of solution methods.

Table 3 :
Results for the undamped cantilevered Euler-Bernoulli beam in Abu-Hilal (2003) with intermediate simple support and subjected to a concentrated sinusoidal force applied at the free end.

Table 4 :
Abu-Hilal (2003)ournal of Solids and Structures 11 (2014) 2171-2202 Results for the damped cantilevered Euler-Bernoulli beam inAbu-Hilal (2003)with elastic support at the free end and subjected to a concentrated sinusoidal moment applied at intermediate location.

Table 5 :
Foda and Albassam (2006)ler-Bernoulli beam inFoda and Albassam (2006)with an intermediate simple support at various locations and subjected to a concentrated sinusoidal force applied at the free end table, the reliability of the present method could be confirmed by comparison with a FEM model with 1000 beam elements.
a Harmonic analysis by direct integration with 1000 BEAM3 elements and consistent mass matrix

Table 6 :
Undamped simply-supported Timoshenko beam in Foda and Albassam (2006)] Subjected to a concentrated sinusoidal force applied at half span and with (a) attached mass and (b) two symmetrically attached springs Latin American Journal of Solids and Structures 11 (2014) 2171-2202