Flutter Instability of Timoshenko Cantilever Beam Carrying Concentrated Mass on Various Locations

This paper presents effects of shear deformation on flutter instability of cantilever beam subject to a concentrated follower force. The discrete form of equation of motion is formulated based on the Lagrange. In the presented formulation, the beam is modeled using Timoshenko beam theory, and constant shear distribution through the thickness of the beam is considered. Consistent interpolation scheme is adopted to avoid the shear locking for thin beams. Consequently, convergence of the finite element simulation is enhanced. The effect of rotary inertial term is considered in the flutter study, which has significant influence on the beam behavior as the beam thickness increases. The axial degrees of freedom are taken into account in energy expressions, to improve the accuracy of the results. Results presented for different beam geometries. The numerical results show high efficiency and good convergence characteristic. The effect of concentrated mass on the flutter instability of beam is considered and results are presented for various locations and values of concentrated masses. Furthermore, the shear effects are highlighted in this study by comparing the results obtained from the Euler-Bernoulli with those obtained from the Timoshenko beam model.


Latin American Journal of Solids and
beam.Forces are conservative if their components are the negative derivatives of a scalar potential function, and dependent of spatial coordinate only, such as gravitational force, which is a wellknown sample of conservative forces.In some cases, such as those studied in this paper, some parameters of force acting on the system are a function of system parameters.In this case, the force is dependent on the system parameters and is no longer conservative.The clamped beam under the action of tip follower force is known as the mathematical model for special aerospace structures (Sugiyama, 1999 andMallik et al. 2005).In addition, as a suitable model for investigating behaviors of pipes conveying fluid has been considered by many researchers.With some modification, this model is applicable to the vibration problems of automotive disk and drum-brake systems that involve dry friction (Mottershead et al., 1997 and1995).For systems consisting of non-conservative forces, exact solutions (Beck, 1952 andBolotin, V. V. 1965) are complicated or do not exist.Therefore, stability analysis regarding these systems is mostly based on approximate techniques, such as finite element method, Galerkin method, etc. Extensive survey of the problem involving nonconservative forces is presented in Elishakoff (2005) and Langthjem et al. (2000) where the history of analysis and controversies about the Beck's column are presented in a comprehensive manner.The dynamic analysis of cantilever beam under follower force and intermediate spring support is investigated by Lee (1995) based on the Lagrange approach solved using the assumed modes.In his paper, Lee employed the Euler beam theory.The applied mode shapes are the same as vibration mode shapes of a cantilever beam, obtained from solution of associated differential equation along with boundary conditions (Clough and Penzien, 1975).He also extended the method to the tapered cantilever beams on Winkler-type elastic foundation and proposed a method to take in to account the effect of viscous damping of foundation on the flutter instability of cantilever beam (Lee, 1996).There is a considerable difference between the vibration mode shapes of cantilever beam in flutter state and the natural modes of the cantilever beam.Consequently, Lee's solution shows a substantial error in the vibration analysis of beam under a follower force.The Euler beam theory, which assumes zero shear stress through the thickness of the beam, is not applicable for stubby and short beams, where transverse shear stress can severely affect the vibration behavior of beam.The Euler beam theory underestimates the deflection and overestimates the natural frequencies and buckling loads of the thick beams.As a result, in designing thick, beam like structures Timoshenko beam theory leads to a more accurate and acceptable results.Dynamic stability of Timoshenko beams are studied by Ryu and Sugiyama (1994).They also studied the effect of concentrated tip mass on the stability of beams under follower force.Moreover, for investigation of the effects of shear deformation, they introduced a factor that contains parameters such as shear correction factor.The formulations are defined by the out of plane deflection only and axial degrees of freedoms are neglected.Ignoring these degrees of freedom cause errors in predictions of thick beams behavior and neglects the effect of axial modes.By extending the method, Ryu et al. (1998) studied dynamic stability of a vertical cantilever beam under sub-tangential follower force.In their study, the Follower force is assumed to be summation of a follower force and a dead load from a rigid body on the end of a vertical column.They used finite element method and presented results for different geometries and conditions of beam.Kim and et al. (2000) investigated the effect of crack on the dynamic stability of a Free-Free beam subjected to a follower force.Wang (2004) studied the effect of the crack intensity and loca-tion on the buckling or flutter compressive load of a beam with a single crack.Viola et al. (2007) considered the effect of sub-tangential forces on the dynamic stability of a cracked beam.Caddemi et al. (2014) investigated the stability of multi-cracked cantilever Euler beam-column subjected to conservative or nonconservative axial loads.
In this paper, the finite element method in context of Timoshenko beam theory is used to predict dynamic instability behavior of cantilever beam subjected to an end follower force.In current studies, non-linear strain terms are used to derive a geometric stiffness matrix, taking in to account the axial degrees of freedom effects.The linear interpolation functions are adopted for the axial displacements and rotation terms, while the quadratic interpolation functions are employed as trial functions for out of plane deflection of the beam.As a result, a consistent interpolation scheme is achieved for the transverse shear strains, and the shear locking for thin beams is avoided and convergence of the analysis is enhanced.Details of the present method are described in the following sections and the components of finite element procedure are discussed in detail.As stated, Beck's column is considered mostly for analysis of hollow shape structures such as pipes, etc.In the case of complicated cross sections, the shear correction factor and other cross sectional parameters are crucial to achieve the accurate analysis.These effects are investigated in the present study.In some of the applications, Beck's column model is applied for analyzing the structures with variable mass and changes in location of mass center.For modeling such structures, it is possible to assume a concentrated mass on a beam where location and magnitudes of mass changes.The results for flutter instability of beams with concentrated masses are presented and a parameter study on the mass value and its location is performed.The results of dynamic analysis are also compared for various beam slenderness ratios.

THE BEAM CONSTITUTIVE RELATIONS
The Timoshenko beam theory (TBT) presented in this paper is based on the same assumptions as those of the Euler beam theory, except that it no longer assumes that the straight planes normal to the neutral axis remain normal after deformation.Instead, it assumes the normal plane rotates with respect to normal to the neutral axis before deformation as shown in Figure 1.Wang, et. al, 2000).Structures 13 (2016) 3005-3021 The displacement fields for the TBT (u; w) in Cartesian coordinates (x; z) where behavior of the beam after deformation is defined can be expressed as follows (C.M. Wang, et. al, 2000):

Latin American Journal of Solids and
where 0 u , 0 w and x f are the unknown functions which define the displacement relations of Timoshenko beam.It is possible to convert this equation into the Euler beam theory by a simple substitution.The condition ,x x w f = -in equation ( 1), yields the same displacement field as that of the Euler beam theory (EBT).The variational functional of the above displacement field requires the satisfaction of C 0 continuity of variables.The degree of continuity for parameters must be considered in choosing interpolation function for each unknown function in finite element formulation.
Employing the linear strain-displacement relationships of three-dimensional elasticity theory in conjunction with equation (1), the five strain components of TBT at a general point can be expressed in terms of the fundamental neutral axis quantities.Thus, the in-plane direct and shear strains of TBT are defined as follows: From the equation ( 2) where the strains of beam are defined, it can be noticed that normal strain of the beam is a linear function of vertical displacement coordinate Z (same as Euler beam, theory) and the shear strain is defined as a constant value.The shear strain of the beam from the elasticity and equilibrium equation of beam is obtained as a quadratic function of the vertical displacement.To compensate the error caused by assuming the constant variation of shear strain through the thickness, shear correction factor is used in this analysis.
The constitutive equations for the beam, the geometry is shown in Figure 1, can be obtained through the use of equation ( 2), Hook's law and appropriate integration through the thickness of the beam.The constitutive equation which defines the relation between the stress and strain of the beam in a matrix form is expressed as: where Mx, Nx and Qx are the stress resultants of the beam which is defined by integration through the thickness of beam in equation ( 5) and IIi is described in equation ( 6): Latin American Journal of Solids and Structures 13 (2016) 3005-3021 (5) Parameters for inertia moments of cross section which appeared in the constitutive equation of the beam are defined as below: In the above equation, d is the width of beam which is a specific height and defined as a function of vertical coordinate (z).
Shear correction factor, ks, also appeared in the constitutive equation.The calculated shear strain through thickness of the beam by using equilibrium stresses, assumes parabolic distribution of stress through the thickness.In contrast, the Timoshenko beam theory assumes the constant state of shear stress through the beam height.To take into account the difference of energy due to shear, which is predicted by either of these theories, the ks shear correction factor is introduced (C.M. Wang, et. al, 2000).This factor is defined using to material and geometry of the beam.For beams with rectangular and circular cross sections and isotropic material, it is assumed 5/6 and 9/10, respectively.Detailed formulation for shear correction factor of shallow circular beams will be described later in this paper.In constitutive equation of beam, E is the Young module of elasticity and G is termed as shear module.Visual descriptions of parameters defining the geometry of the beam are given in Figure 2.For small deformations of beam, the angle between the deformed and un-deformed neutral axis of the beam configuration at the end is assumed to be identical to the rotation φx.

ENERGY EXPRESSION OF A BEAM ELEMENT
The extended Hamilton's principle for small motions of the beam under tip follower force is expressed as follows: In the above equation, U is the strain energy of the beam element, T is the kinetic energy for vibration motion of the beam, Vis the loss in potential energy or the potential energy due to conservative component of follower force, andWnc is the virtual work, which resulted from the nonconservative component of tip follower force.The strain energy of a beam element is defined with respect to the material properties of the beam and the strain vectors as follows: e and Dm are the strain vector and material matrix of the element respectively.On substituting Equations ( 2) to (3) and then into Eq.(8), it is possible to obtain the strain energy of an element as a function of unknown displacement functions: The kinetic energy of beam in this study is defined as a summation of the kinetic energy due to vibration of the whole beam and the movements of concentrated mass, CM, which is placed in a different arbitrary location of beam as shown in Figure 2. The kinetic energy of a vibrating beam is defined as: In the above equation, ρ, is the material density of a beam which is assumed to be constant in the beam's length.Using equation (1) and integrating through the thickness of the beam, the kinetic energy of beam element is defined as: In equation above TCM is dealing with kinetic energy due to the concentrated mass, which can be a mathematical model for heavy machinery or any abrupt change in the density of beam.This energy is neglected in all elements except one, which contains the concentrated mass.For this element, the expression below is added to the kinetic energy of the beam element: Latin American Journal of Solids and Structures 13 (2016) 3005-3021 In the present formulation and in order to obtain the characteristic matrices, follower force on the beam is divided into two separate components; one with constant magnitude whose direction always remains parallel to the neutral axis of the un-deformed beam.This component is conservative and the potential energy due to this conservative force is defined as: By importing the displacement functions of TBT from equation (1) into the equation above the potential energy VC is expressed as: It is clear that in this formulation the non-destabilizing effect corresponding to rotational degrees of freedom is considered in the loss potential energy for the beam.The direction of the second component is normal to the neutral axis of an un-deformed beam and its magnitude changes corresponding to rotation of the normal line.This component of force is a function of the degrees of freedom of the beam.Thus, it is non-conservative and no potential energy is defined for it.Indeed, the virtual work done by this transverse component can be attained.For small deflections of beam, the vertical component of follower force is defined as a product of its magnitude and the rotation of beam at the end.Finally, the virtual force is defined by equation below, which is function of degrees of freedom at the point on which the load is imposed.

THE FINITE ELEMENT FORMULATION OF THE PROBLEM
The unknown parameters in equation (1) are defined as a time dependent function.It is assumed that the unknown degrees of freedoms are combined from the time dependent function and displacement degrees of freedom as the following equation: Considering the above equation, variation of each of the displacement functions is represented spatially as a product of Lagrangian polynomial and degrees of freedom vector. where: Latin American Journal of Solids and Structures 13 (2016) 3005-3021 M and N are linear and quadratic Lagrange polynomials, respectively.Visual description of these functions is shown in Figure 3. Consequently, for the purpose of interpolation, two nodes are needed for linear shape function and three nodes for quadratic function.As a result, the beam elements will have three node and seven degrees of freedom each.In Figure 4, b corresponds to the length of an element and described as: where, ne is the number of elements for the beam analysis.Approximation for unknown functions in equation ( 18) causes φx and w,x to be the polynomials from the same degree.Consequently, the consistent interpolation (N.Reddy, 2006) is applied in the present study.As a result, shear locking is avoided and convergence of method is enhanced, especially for the analysis of thin beams.
By substituting the element displacement fields given in equation ( 17) into equations ( 8)-( 15), expressions can be derived for the strain energy, kinetic energy, potential energy of the applied load, and also virtual work of the non-conservative component of follower force.These expressions are quadratic functions of the complete set of displacement degrees of freedom of each element.For each element, the column matrix of {d} for total degrees of freedom of beam element is defined as below: Finally, the expressions for energy and work will be written as the familiar forms of: In the present study, these expressions are obtained using differentiation with respect to each degree of freedom.The matrices are containing polynomials as their elements.Thus, it is possible to carry out integrations involved in evaluating the matrices exactly over the length of the element, which makes the solution of the problem faster.In addition, it is possible to use methods such as gauss integration in finite element formulation.
Obtaining natural frequencies of structures generally with respect to imposed end force requires one to solve an eigenvalue problem as follows: where  ˆ is the corresponding eigenvector and e K , G K and a M are the elastic stiffness matrix, geometric stiffness matrix and mass consistent matrix of the whole beam structure respectively.
In the present study, G K is a nonsymmetrical matrix due to non-conservative force and is defined as a combination of a symmetric geometric matrix (due to conservative component of follower force) and the non-symmetric matrix according to its non-conservative component as below: These stiffness matrices can be assembled in the standard fashion by defining a simple connectivity matrix (N.Reddy, 2006).Equation ( 23) is solved using an appropriate eigenvalue solution procedure that will give the beam's natural vibration mode shapes and frequencies at the applied follower force (Bathe, K. J. 1982).Structures 13 (2016) 3005-3021 5 NUMERICAL STUDY Ryu et al. (1994) investigated the effect of shear deformation on flutter stability of the beam through introducing the shear deformation parameters which are defined as a function of beam geometry as below:

Latin American Journal of Solids and
To validate the obtained result from the current method they are compared with those of Ryu et al. (1994), for different values of shear deformation parameters and the results are shown in Table 1.The results are obtained for rectangular cross section with ks=5/6 for different shear deformation factors and corresponding slenderness ratios.It should be noted that in a comparison between results there exists a small difference in the results from the present method and those reported in Ryuet al. (1994).Destabilizing effect of shear degrees of freedom (equations ( 13)-( 14)) was not considered in Ryuet al. (1994),and analysis was based only upon the lateral deflections.It is clear that neglecting these degrees of freedoms in geometric matrix of a beam results in a slight error of the analysis, especially for stubby beams.The values of critical flutter load for different slenderness ratios are presented in Table 2 Analysis is based upon both the Euler-Bernoulli and the Timoshenko beam theories.Results are compared with those reported by Lee (Lee, 1996) obtained using assumed modes method and the exact solution of differential equation of Euler-Bernoulli beam under follower force, by applying dynamic instability criteria (Beck et. al 1952, Bolotin, V. V., 1965).It is worth to note that difference percentage of these result are presented in bracket in Table 2.As expected from these analysis, by increasing the beam thickness, the effect of shear deformation increases and there is more difference between the results from the Euler formulation and the present formulation, which accounts for shear deformations.As one may notice for beams with higher slenderness, the results are identical to those of Euler beam theory, as expected.
The effect of slenderness ratio on variation of natural frequencies due to follower force is shown in Figure 5 where the changes in first and second natural frequencies of cantilever beam are shown for thin and thick beams.The non dimensional parameters are defined as Pb = (PEII2)/L 2 and Wb=ω(ρAL 4 /EII2) 1/2 .By increasing the follower force to the specific value, the paths for frequency alteration due to the follower force coalesce and the beam flutter occurs.Effect of follower force on first and second mode shapes of a cantilever beam is shown in Figure 6.It is clear that by increasing the magnitude of follower force, mode shapes are merging to each other and reach or coalesce into one mode shape after critical follower force.
As seen from equation ( 11), (kinetic energy of beam element), the presence of rotary inertia effect cause to increase the kinetic energy of system, as expected.It is clear that the natural frequencies of beam decreases with increasing of kinetic energy.Consequently, reducing natural frequencies of system causes flutter boundary is moved left, i.e., the beam flutter occurs quickly.Figure 5 and results of reference (Park, 1987) show the correctness of the above content.
Beams with circular cross sections or pipe shape structures are more popular in various industries.As a result studying of these types of cross sections is of paramount importance in predicting behavior of pipes and aerospace structures (Curtis Petrus, 2006).The effect of a cross section on flutter stability is considered in Table 3, where values of a critical flutter load are presented with respect to the changes in the cross section type.Shear correction factor in circular cross section is a function of inner and outer radius of hollow beam (Cowper, 1966): where S is the radius ratio and defined as ratio of inner Ri radius to outer radius Ro of hollow beam: The change in the cross section leads to the variation of shear correction factor, which is defined as a function of radius and Poisson ratio in equation ( 25).As noted before, the effect of shear deformation effect is not negligible in analysis of thick beams.As a result, the changes in the values of critical follower force due to change in cross section (in Table 3) is more apparent for thick beams.
Effect of tip mass on flutter stability is presented in Table 4.The effect of tip force on flutter stability in the present method is compared with Ryu et al. (1994) and a good agreement between the results is achieved.The results are presented for different values of mass ratio M*which is defined as ratio of tip mass to the mass of beam.This study is for analyzing a beam response with concentrated mass (CM), in which xCM is equal to the length of the beam.The effect of position of concentrated mass on natural frequencies of cantilever beam is depicted in Figure 7 where results from the present method are compared with the results from Dunkerley's Latin American Journal of Solids and Structures 13 (2016) 3005-3021 approximation technique (Tse et al., 1978) and good agreement can be observed.It is clear from the results that the natural frequency of a cantilever beam decreases with increase of μ.
Figure 7: Effect of μ on first natural frequency of cantilever beam with concentrated mass (μ=xCM/L).
Figure 8 shows the critical follower force of a cantilever beam (Flutter point) with concentrated mass for different values and location of mass.It is clear from the results that the concentrated mass on the middle of cantilever beam increases the flutter force of a beam for all shear deformation parameters.Moreover, concentrated mass on the first or last quarter of a beam has a destabilizing effect on the beam in case of flutter.

CONCLUSIONS
In this paper, the finite element method is applied to flutter analysis of a cantilever Timoshenko beam under non-conservative follower force.The effect of shear deformation and rotational kinetic energy are taken into account.Results are compared with those reported elsewhere, and comparison among the results is presented.The effect of concentrated mass on flutter instability of cantilever mass is inspected and has been shown that the concentrated mass increases or decreases the critical flutter force of beam depending to its location on the beam.The effect of axial and rotational degrees of freedom inertias in Timoshenko beam is considered in the present study.It is concluded from the results that including the inertia of these degrees of freedom causes a changes in the result for critical follower force.In addition, it can be concluded from the obtained results that the effect of shear deformation in the critical follower force cannot be disregarded specially in the case of stubby beam cases.

Figure 2 :
Figure 2: Geometry of cantilever beam under follower force

Figure 5 :
Figure 5: Variation of first and second natural frequency with respect to applied follower force for different slenderness ratios.

Figure 6 :
Figure 6: Mode shapes of a cantilever beam under follower force for different values of follower force.

Figure 8 :
Figure 8: Critical follower force of cantilever beam for different values and locations of concentrated mass.

Table 1 :
Effect of shear deformation parameters on critical follower force Pn = Pcr E II2/L 2

Table 2 :
. Effect of slenderness ratio on critical follower force.

Table 3 :
Effect of radius ratio on critical follower force of stubby and slender beams.