Dynamic Instability of Beams Under Tip Follower Forces Using Geometrically Exact , Fully Intrinsic Equations

In this study, the dynamic instability of beams under tip follower forces are considered. The beam is modeled by using the geometrically exact, fully intrinsic beam equations which is subjected to an inclined tip follower force. Generalized differential quadrature method is employed to solve the governing equations. The effect of different parameters such as follower force inclination and magnitude, rotating speed, the distance between the beam center of gravity and elastic center, and cross-sectional properties on the instability boundary of beams are examined. Numerical results reveal that the critical load of the system can be influenced or in some cases be reversed by the combination of these parameters instead of considering these parameters separately. Moreover, it is shown that notwithstanding the simplicity of the equations and the method of solution, the results are very accurate and therefore the fully intrinsic equations and the method of solution is very useful for the dynamic solution of rotating and non-rotating beams.

Latin American Journal of Solids and Structures 13 (2016) 3022-3038 structures exposed to follower forces, as mentioned by Bolotin (1963), was the work done by Nikolai, and after that many papers and also some books are concentrated on this topic (Bolotin (1963(Bolotin ( , 1964)), Leipholz (1978), Simitses and Hodges (2006)).In the spite of all the published works, there seems that very little attention has been put on the lateral-torsional stability of rotating cantilever beams under inclined transverse follower forces.In 1952, Beck (1952) determined the dynamic instability of a cantilever column exposed to a tangential follower force.In this study, the follower force was a concentrated force that exactly applied at the tip of the column.Leipholz (1962) considered the stability of a cantilever column under uniformly distrusted tangential follower force approximately.Como (1971) studied the stability of a cantilever beam subjected to a lateral tip follower force without considering the mass and inertia distribution of the beam.In this study, a concentrated mass and inertia located at the tip of the beam were included.This work has been reconsidered by Wohlhart (1971) by considering the mass and inertia distribution of the beam, and the effect of different parameters on the stability of the beam was presented.Zuo and Schreyer (1996) studied the dynamic and static instability of cantilevered beams as well as simply supported plates subjected to non-conservative forces.The lateral stability of a beam subjected to follower forces has been investigated by Detinko (2002)).In this paper, it was shown that by neglecting the slight internal and realistic external damping, the critical load of the system will be determined inaccurately.Feldt and Hermann (1974) considered the bending torsional instability of cantilever beams under tip transverse follower forces.They presented a comprehensive study on different parameters on the instability boundary of the system but the obtained results were not in agree with previously published works.The effect of shear deformation and rotary inertia on the stability of slender Euler-Bernoulli cantilever columns subjected to follower forces has been considered by Nair et al. (2002).Stability determination of a rotating cantilever subjected to dissipative, aerodynamic, and transverse follower forces has been investigated by Anderson (1975).In this study, the variation of critical flutter load with respect to the hub radius, the rotational speed, aerodynamic load parameters and the warping rigidity has been examined, but the effect of follower force inclination angle and the mass centroid offset from the reference line has been not determined.Hodges (2001) investigated the lateral-torsional flutter of a deep cantilever beam subjected to a tip lateral follower force.In this paper, the effect of different parameters on the critical follower force has been studied.This work has been continued by adding the aerodynamic loading on the beam (Hodges and Patil (2002)).More recently, Fazelzadeh and Kazemi-Lari (2014) and Kazemi-Lari and Fazelzadeh (2015) considered the stability analysis of a deep cantilever beam with transverse uniformly and partially distributed follower forces.They observed that the position and magnitude of the distributed load can influence the stability boundary of the beam.The effect of engine thrust that may be modeled by a transverse follower force on the aeroelastic instability of wings has been considered by Hodges and Patil (2002), Amoozgar et al. (2013), Fazelzadeh et al. (2009), Mardanpour et al. (2013) and Mardanpour et al. (2014).In all these studies, it was shown that the magnitude and position of the follower force influence the stability boundary of the wing.
In all the studies mentioned above, the beam has been modeled by displacement-based beam theories and to the best knowledge of the authors, no one used the stress based geometrically exact beam theories to study the stability boundary of rotating beams under follower forces which is very accurate.On the contrary of the displacement-based theories, some theories are developed in which Latin American Journal of Solids and Structures 13 (2016) 3022-3038 neither displacement nor rotation variables appear in the related partial differential equations which are stress (or strain) based formulation.These theories are called fully intrinsic theories based on the terminology of Green and Laws (1966) and Reissner (1973) and originally developed by Hegemier and Nair (1977) and more recently reconsidered by Hodges (2003).In recent years, many researchers considered this type of beam equation to study the static and dynamic behavior of beams.Patil and Hodges (2006) studied the flight dynamic of a flexible flying wing configuration by using the fully intrinsic geometrically exact beam equations.Mardanpour and Hodges (2013) examined the effect of engine on trim and aeroelastic stability of flying wing vehicles.They showed that considering the gravity can influence the stability envelope of the wing.This work has been continued by considering multiple engines and time-dependent thrust effects on the aeroelastic behavior of flying wings.To the best of the author's knowledge, the dynamic stability of rotating cantilever wings subjected to inclined follower forces by using the geometrically exact fully intrinsic beam equations have never been considered and the objective of the present paper is to consider this topic and further investigation of different parameters on the stability boundary of these beams.
Up to now, there are a number of solution methods that are used for solving the fully intrinsic beam equation.They are limited to simple discretization of the equations (Sotoudeh et al. (2010) and Mardanpour et al. (2014)), variable order finite element method (Patil and Hodges (2011)), Galerkin method (Patil and Althoff (2011)), and Chebyshev collocation method (Khaneh Masjedi andOvesy (2014, 2015)).All these methods are successfully applied to the fully intrinsic beam equations for different problems.In this paper, an alternative method of solution called Generalized Differential Quadrature (GDQ) method will be used to solve the fully intrinsic beam equations.This method is among the easiest and efficient solution methods which can be applied directly to the partial differential equations.The differential quadrature method was first introduced by Bellman and Casti (1971).In this method, the derivative of a function at a specific point is approximated as a weighted linear summation of the values of the function at all other sampling points along the domain.Shu and Richards (1992) proposed a method to overcome the DQ drawbacks by introducing the GDQ method based on the analysis of a polynomial vector space.By using the GDQ method, the first order derivatives are computed by a simple algebraic formulation without any limitation on choosing the grid points.Until now, many researchers used this method for solving various problems containing partial differential equations in fluid flow (Shu and Richards (1992), Shu et al. (1995) and Shu et al. (1996)) or in structures (Bert and Malik (1996), Du et al. (1995), Du et al. (1996), Laura andGutierrez (1993, 1994) and Lin et al. (1994)) and the literature is overwhelmingly rich, and here only a few examples are presented.Lin et al. (1994) used the GDQ method to study the plate deflection with nonlinear boundary supports.Marzani et al. (2008) solved the nonconservative stability problems via GDQ method.Lal and Saini (2015) determined the vibration behavior of non-homogeneous orthotropic rectangular plates of bilinearly varying thickness by GDQ method.Also, a review on differential quadrature method and its applications in computational mechanics was implemented by Bert and Malik (1996).Amoozgar and Shahverdi (2016) used the generalized differential quadrature method for solving the geometrically exact fully intrinsic beam equations and showed that this numerical method is an efficient method for these type of equations.
In this paper, the instability of beams under tip follower forces modeled by geometrically exact fully intrinsic beam equations is investigated by using the generalized differential quadrature meth-Latin American Journal of Solids and Structures 13 (2016) 3022-3038 od.First, the governing equations of the system with tip follower force have been demonstrated.Then by introducing the method of solution to the governing equations, the beam bending and torsional modes are obtained.The effect of rotating speed, the inclination angle of the follower force, the offset of the mass centroid from the elastic axis, the ratio of the bending and torsional frequencies and some other parameters on the stability boundary of the wing are examined.It is found that the combination of these parameters can change the instability region of the system dramatically.

GOVERNING EQUATIONS
Figure 1 shows a beam with its reference coordinates in deformed and undeformed configurations which are denoted by b(x1), and B(x1,t), respectively.The nonlinear geometrically exact, fully intrinsic governing equations for the dynamics of an initially curved and twisted beam undergoing large deformations in the deformed coordinate can be expressed as (Hodges (2003)): where ( ) is the partial derivative with respect to x1, ( )  shows the partial derivative with respect to time.M and F are the internal moment and force vectors which are expressed in the matrix form.H and P are the generalized angular momentum and linear momentum vectors, Ω and V are the generalized linear and angular velocity f and m denote the generalized strains and curva- ture, κ and γ are the external moment vectors and force applied on the beam, respectively.It is noted that all above variables are in the matrix form and expressed in the deformed state except the initial curvature k wich is expressed in the undeformed configuration e1 is a vector wich its arrays are as bellow: The generalized strains of an isotropic beam are related to the generalized force and moment via the following constitutive equations as follows (Simitses and Hodges (2006)): On the other hand, the generalized linear and angular velocities are related to the generalized linear and angular momentum through the cross-sectional inertia matrix as: where , I  and  are the cross-sectional inertia matrix, the mass per unit length, and the mass center offset from the beam reference axis, respectively.The cross-sectional inertia matrix and the mass center offset vector components are: where i2 and i3 are the cross-sectional mass moments of inertia and i23 denotes the cross-sectional product of mass moment of inertia, and 2 x and 3 x are offsets from the reference line along the cross-section coordinates.By inserting Eqs. 3 and 4 in the fully intrinsic equations of motion (Eq.1), the resultant equations can be rearranged in a manner that the primary unknown variables will be the generalize force, moment, linear velocity and angular velocity.There will be 12 unknown variables and therefore 12 boundary conditions are needed for solving the fully intrinsic equations of motion.In this paper, a cantilever beam with length L is considered and the 12 boundary conditions will be: It is of note that, the displacements and rotations variables of the beam can be determined by the following relation in post-processing phase (Sotoudeh et al. (2010)): where u is the displacement vector of the beam.

NUMERICAL SOLUTION PROCEDURE
The generalized differential quadrature method is used to solve the fully intrinsic equations.The main idea of this method is to evaluate the derivative of a function at a specific point as a weighted linear summation of all the functional values at all other points along the domain (Shu and Richards (1992)).In this method, the n-th order partial derivative of a function like f(x) with respect to the space variable can be written as: where N is the number of sampling points considered in the domain and W is the weighting coefficients which can be written as follows (Shu and Richards (1992)): Latin American Journal of Solids and Structures 13 (2016) 3022-3038 (1) (1) (1) ( ) , , 1, 2,...., N and (x x ) M ( ) where Moreover, the weighting coefficients of the second and higher order derivatives will be: ( 1) ( ) ( 1) (1) , 1, 2,....., 2,3,...., 1 , 1, 2,... , 1, 2,..., 1 It is noted that the above weighting coefficients depend on the derivative order and the number and distribution of sampling points.In this study, the Chebyshev-Gauss-Lobatto point distribution as shown in Figure 2 is used (Shu (2000)): ( 1) 1 cos 2 (N 1 ) where is the length of the domain (beam).By using the GDQ method, the primary variables of the equations of motion can be discretized as: (1)  and Structures 13 (2016) 3022-3038 (1) By substituting Eq. 16 into the equations of motion (Eq.1), the discretized form of the differential equations will be: Furthermore, the discretized form of the boundary conditions for the clamped beam is: .
By rearranging the equations of motion, the resultant equations in the matrix form can be written as: Where, A and B are the matrixes of linear coefficients, C is the matrix of nonlinear coefficients, D is the vector of external forces and moments, and q is the vector of unknown parameters.It is of note that the total system of equations consists of 12N equations and unknowns.The Eq. 19 is a set of nonlinear equations and the solution of the system consists of two steps.First, the nonlinear steady-state solution of the equations must be determined, and in the second step, the equations must be linearized about this nonlinear steady-state solution.By dropping the time derivatives terms of the system the nonlinear steady-state solution of the system can be obtained by using the Newton-Raphson iterations.The resulting equations can be written as: where q is the steady-state solution because of the steady-state forcing determined in D .By using this steady-state solution, the linearized equations of the system can be obtained which can be used to determine the eigenvalues of the system.The linearized equations may be written as: Latin American Journal of Solids and Structures 13 (2016) 3022-3038 To form the above set of equations into a generalized eigenvalue problem, it is assumed that: ˆexp( ) where is the eigenvalues of the system.By using the above equation, the resulting equation is in the form of a standard eigenvalue problem as: By using the Eq. 23, the eigenvalues of the linearized system can be determined.

NUMERICAL RESULTS
In order to check the validity of the developed code, the obtained results are compared with those reported in the literature.It is worth mentioning that the authors have been checked the efficiency of the proposed GDQ method of solution with respect to the conventional FEM method in their previous paper (Amoozgar and Shahverdi (2016)).For numerical usefulness, the following nondimensional parameters are introduced: where  is the cross-sectional mass radius of gyration, S is eigenvalues of the system, e is the mass centroid offset from the reference line in the B2 direction, and 1  and 1  are same as those used in Hodges (2001).The variation of the imaginary parts of the eigenvalues of the system, versus K for e=0 and r=3.8 and r=2/3 are calculated and compared with those reported in literature in Figs. 3 and 4. It is of note that in these two plots due to the independence of the eigenvalues to σ for e=0 (Hodges (2001)), the value of σ is not mentioned.
As it is clear from these figures, the obtained results are in good agreement with the results reported by Hodges (2001).On the other hand the variation of the critical force for two different values of e and for 0.05   is illustrated in Figure 5.By taking a quick look, it is realized that the obtained results have a very good correlation with those presented in the literature.
In the following, it is considered that the tip follower force as shown in Figure 6    , versus the inclination angle is plotted.In this case, by increasing the inclination angle, the critical force of the beam increases.This trend is about the same for all values of e.Notice that by setting the inclination angle to / 2    , the critical force of the system reaches to the Beck's column critical force (Beck (1952)) which for e=0 is obtained as Kcrit=406.In a similar manner of Figure 7, the variation of the critical force versus inclination angle for three different values of r and for e=0.03 is depicted in Figure 8.It is observed that for all values of r, by increasing the inclination angle, the critical force increases and the manner of increase is the same for all values of r.To examine the effect of rotating speed on the critical force of the system, the beam is considered to rotate about its third axis(B3) and simultaneously subjected to a tip follower force.Figure 9 shows the variation of the critical force against the rotating speed of the beam for various values of r and for e=0.03 and 0.05

 
. As the rotating speed increases for r=1, the critical load of the beam decreases monotonically until reaching to zero, but the behavior of the system for other two values of r is a little different.One finds that for r=1.5 r=2, by increasing the rotating speed, the critical load decreases, but the rate of reduction is variable.For example for r=2, first, the critical load reduces slowly until 3   , and from this point until 8   , the critical force decreases rapidly and after that again it decreases slowly.The variation of the critical load by increasing the angular speed of the beam for two different values of e and for r=1.5 is calculated and plotted in Figure 10.It is concluded that the manner of variation of critical load versus rotating speed is surprisingly different for e=0 and e=0.03.In the other words, for e=0, by increasing the angular speed, the critical load increases while for e=0.03 this trend is inverse.It is noted that the results for e<0 is not presented here because for larger values of rotating speed, the imaginary parts of the eigenvalues do not coalescence to each other and they simply veer each other as shown in Figure 11 for e=-0.03 and Ω=30.
The effect of inclination angle on instability boundary of rotating beams for various values of rotating speed is demonstrated in Figure 12 for e=0, r=1.5 and 0.05

 
. It is found that for all values of rotating speed, by increasing the inclination angle, the critical load increases and this trend is almost the same for all values of rotating speed.

CONCLUSION
In this paper, the instability of a rotating and non-rotating beam under an inclined tip follower force has been considered.The beam is modeled by using the geometrically exact fully intrinsic beam equations.The generalized differential quadrature method is used to convert the partial differential equations to ordinary ones.The stability of the system is determined by investigation the eigenvalues of the linearized system about its equilibrium state.It was found that the proposed method of solution simulates the system accurately.Moreover, the variation of the critical load versus different parameters has been examined, and based on the numerical results, the critical load of the system can be influenced by changing the values of the follower force inclination and magnitude, rotating speed, distance between the beam center of gravity and elastic center, and crosssectional properties.Finally, it can be highlighted that the geometrically exact fully intrinsic beam equations are very useful for dynamic instability of rotating and non-rotating beams subjected to tip follower forces.

0R
and r are the reference line position vectors in the deformed and undeformed states, re- spectively and u is the displacement vector.

Figure 3 :
Figure 3: Variation of the imaginary parts of the eigenvalue versus K for r=3.8 and e=0.

Figure 4 :Figure 5 :
Figure 4: Variation of the imaginary parts of the eigenvalue versus K for r=2/3 and e=0.

Figure 6 :
Figure 6: Schematic of the cantilever beam subjected to an inclined follower force.

Figure 7 :
Figure 7: Variation of critical load versus inclination angle for r=1.5 and

Figure 8 :
Figure 8: Variation of critical load versus inclination angle for e=0.03 and

Figure 9 :
Figure 9: Variation of critical load versus nondimensional rotating speed for e=0.03 and

Figure 10 :
Figure 10: Variation of critical load versus nondimensional rotating speed for r=1.5 and