On the Flexural-Torsional Vibration and Stability of Beams Subjected to Axial Load and End Moment

The free vibration of beams, subjected to a constant axial load and end moment and various boundary conditions, is examined. Based on the Euler-Bernoulli bending and St. Venant torsion beam theories, the differential equations governing coupled flexural-torsional vibrations and stability of a uniform, slender, isotropic, homogeneous, and linearly elastic beam, undergoing linear harmonic vibration, are first reviewed.The existing formulations are then briefly discussed and a conventional finite element method (FEM) is developed. Exploiting the MATLAB-based code, the resulting linear Eigenvalue problem is then solved to determine the Eigensolutions (i.e., natural frequencies and modes) of illustrative examples, exhibiting geometric bending-torsion coupling. Various classical boundary conditions are considered and the FEM frequency results are validated against those obtained from a commercial software (ANSYS) and the data available in the literature. Tensile axial force is found to increase natural frequencies, indicating beam stiffening. However, when a force and an endmoment are acting in combination, themoment reduces the stiffness of the beam and the stiffness of the beam is found to be more sensitive to the changes in the magnitude of the axial force compared to the moment. A buckling analysis of the beam is also carried out to determine the critical buckling end moment and axial compressive force.


Introduction
Beams are important and versatile structural elements as many civil, mechanical, and aerospace structures are commonly modeled as preloaded beams or beam assemblies [1].In the past, many researchers have investigated the vibration and stability of various beam structures subjected to diverse loading and boundary conditions.As a result, the governing differential equations of motion for various beam configurations have been already developed.Therefore, from past literatures, it is known that the tensile and compressive axial loads affect the flexural stiffness of beams and increase or decrease the uncoupled bending frequencies and critical buckling loads of such elements.Apart from the above, various coupled vibrational behaviors of beam structures, caused by different geometric [2][3][4][5][6][7][8][9][10] and material couplings [11][12][13][14][15][16][17][18][19][20], have also been thoroughly studied and their associated findings have been well established.The coupled vibrations and stability of axially loaded, or centrifugally stiffened, beam elements have also been investigated [21][22][23][24][25][26][27][28][29][30][31].However, the literature on the effects of combined axial load and end moment and the couplings resulting from the end moment on the stability and vibration characteristics of beams are scarce [32][33][34][35][36][37][38].
Hashemi and Richard [2] used the frequency dependent dynamic finite element (DFE) method, a hybrid method developed in [3] that combines the accuracy of analytical methods to the versatility of numerical methods, to conduct a vibration analysis of beams that are geometrically coupled in bending and torsion.An exact method has been used to determine the flexural-torsional vibration characteristics of a uniform beam with single cross-sectional symmetry by Dokumaci [4].The classical finite element method (FEM) was used by Mei [5] to study the coupled vibration of thin walled beams with an open section.This study included the effects of warping stiffness.The flexural-torsional vibration of a uniform beam was studied by Tanaka and Bercin [6] by determining the exact solution of the governing differential equations.An analytical dynamic stiffness matrix (DSM) solution was formulated by Banerjee [7] for coupled bendingtorsion beam elements.Banerjee et al. [8] also developed an exact DSM to investigate the vibration of an Euler-Bernoulli thin walled beam and exploited the Wittrick and Williams [9] root finding algorithm to extract the Eigensolutions.Once again, Banerjee and Su [10] used the DSM to conduct a free transverse and lateral vibration analysis of a beam coupled with torsion.
Hashemi and Roach [11] also formulated a DFE solution for the free vibration of an extension-torsion coupled composite beam.A quasi-exact DFE formulation for the free vibration analysis of a three layered sandwich beam, consisting of a thick, soft, low strength and density core and two-face layers made of high strength material, was developed by Hashemi and Adique [12].Borneman and Hashemi [13] developed a DFE for the free vibration analysis of bending-torsion coupled laminated composite wing beams.Bannerjee and his coworkers have used the frequency-dependent DSM method for the vibration analysis of isotropic [14], sandwich [15][16][17], and composite [18] beams.Borneman et al. [19] also used the DSM method to investigate the vibrational characteristics of a doubly coupled (material and geometric) defective composite beam.Furthermore, a DSM formulation was presented by Hallauer and Liu [20] to determine the vibrational characteristics and generalized masses of an aircraft wing modelled as a series of three simple beams.
Hashemi and his coauthors [21] developed a DFE formulation to analyse the free vibration of centrifugally stiffened (rotating) beams.Furthermore, Hashemi and Richard [22] formed a DFE solution for the free vibration analysis of axially loaded bending-torsion coupled beams.An axially loaded isotropic Timoshenko beam coupled in bending and torsion was studied by Banerjee and Williams [23].Leung [24] developed an exact DSM of a thin walled beam.Once again, an analytical solution was formulated by Banerjee and Fisher [25] to model a uniform, axially loaded, cantilevered beam with flexural-torsional coupling as a result of noncoincident shear and mass centers.The effects of warping have been neglected in this study.Jun et al. [26] examined the coupled flexural-torsional vibration of an axially loaded thin walled beam with monosymmetrical cross sections by including effects of warping.The effect of axial load has also been previously studied by Murthy and Neogy [27] for clamped and pinned boundary conditions, as well as by Gellert and Gluck [28] for cantilevered boundary conditions.Bokaian [29] determined the natural frequencies of a uniform single span beam subjected to a constant tensile axial load for various boundary conditions.The same author also investigated the vibrational characteristics of a uniform single span beam for ten different end conditions when a constant compressive axial load is applied [30].Shaker [31] conducted a modal analysis to determine the effect of axial load on the mode shapes and natural frequencies of beams.
It has been established by Chen and Atsuta [39] that transverse bending and torsion is coupled by static end moments and that flexural-torsional buckling is comprised of this transverse flexure and axial torsion.Analytical investigations on the influence of axial loads and end moments on the vibration of beams have been previously reported by Joshi and Suryanarayan for a simply supported case [32] and other boundary conditions [33].Joshi and Suryanarayan [34] also studied analytically the flexural-torsional instability of thin walled beams subjected to axial loads and end moments.The same authors examined the coupled bending-torsion vibration of a deep rectangular beam that is initially stressed as a result of the application of moments that vary along the span [35].Pavlović and Kozić [36] developed a closed form analytical solution to investigate the effects of end moments on a simply supported thin walled beam.Furthermore, Pavlović et al. [37] also formulated the analytical solutions to study a simply supported thin walled beam subjected to the combined action of an axial force and end moment.
The reliability and accuracy of such modal analysis results depend on the method implemented.A variety of approaches for the analysis of these structures have been proposed and implemented in the past, including analytical, semianalytical, and numerical methods.Amongst these, the classical FEM method, where beam element matrices are evaluated from assumed fixed (namely, polynomial) shape functions, has been widely used by investigators.This practice results in approximate equations in the form of mass and static stiffness matrices.A deviation from the conventional FEM formulation would pay dividends if improved accuracy of results can be obtained by using shape functions other than polynomials.This is the case when the homogeneous solution of the pertinent differential equation is available for the development of each of the element matrices.For static analysis, the use of the homogeneous solution of the differential equation yields the exact stiffness matrix and load vector for a beam element [40].The use of other alternative shape functions in dynamic FEM formulations has also been explored.
The dynamic stiffness matrix (DSM) method offers a better alternative particularly when higher frequencies and better accuracies of results are required.It relies on a single frequency-dependent matrix which has both mass and stiffness properties.The use of a DSM in vibration analysis is well established [14][15][16][17][18]. Obviously, the method gives more accurate results because it exploits the exact member theory.The matrix is obtained by directly solving the governing differential equation.Other methods have also been explored and reported which are more or less similar to those explained previously.Some of them were developed to treat a particular configuration and group of mechanical systems.
A thorough investigation of the existing conventional FEM and alternative DSM methods in beam vibrations has led to the development of the DFE approach [3].DFE bridges the gap between the standard FEM and the exact DSM methods by advantageously exploiting the generality of FEM and the very precise frequency calculations provided by DSM approach.From this point of view, the method retains the physical aspect of analytical or semianalytical approaches and the power of a numerical method.It has been shown that DFE is an efficient tool for handling periodical structures or systems composed of several identical substructures [3].

Shock and Vibration 3
All methods stated in the above-mentioned references have their inherent advantages and shortcomings.There exists a class of problems for which an exact solution can be obtained.Nevertheless, in most cases, finding a general exact solution for the normal modes and frequencies of the system would be cumbersome, if not impossible, and would involve complicated mathematical procedures.Thus, recourse would be made to one of approximate solutions such as the Rayleigh-Ritz method [41,42] or Galerkin's method [40].
The conventional finite element method (FEM), which uses the Galerkin method of weighted residuals [40], is widely accepted and used for structural analysis.The FEM is adaptable to many complex systems, including those with material and geometric variations, for example, nonuniform geometry.However, it seems that a FEM-based thorough investigation of the geometrically coupled flexural-torsional vibration of beams, simultaneously subjected to both axial force and end moment, and the coupling effects caused by end moments have not been reported in the open literature.Therefore, in what follows, a classical FEM solution of the above problem is presented and the static buckling stability and the coupled free vibrations of preloaded thin (Euler-Bernoulli) beams, subjected to various classical end conditions (with ends free to warp), are investigated.The axial load and end moment are varied and their effects on the beam stiffness and natural frequencies are examined.It is worth noting that, in addition to the conventional beam stiffness and mass matrices [38,40], the presented formulation is characterized by a geometric stiffness matrix, which includes uncoupled and coupled components.The uncoupled component, in turn, is formed by two parts, one associated with the flexure of axially loaded beams [38,43,44] and another relative to the torsion.However, to the best of authors' knowledge, the coupling geometric stiffness matrix resulting from the end moment has not been reported in the literature.The presented FEM formulation is applicable to the members composed of closed sections (e.g., rectangular or square hollow sections), where the torsional rigidity () is very large compared with the warping rigidity (Γ), with ends free to warp, that is, state of uniform torsion, where the twist rate is constant along the span.The presented FEM formulation, however, can also be extended to thin walled beams with open cross sections, where torsion-related warping effects cannot be neglected.

Theory
Consider a linearly elastic, homogeneous, isotropic beam subjected to an end moment, , and an axial load, , undergoing linear vibrations.Euler-Bernoulli bending and St. Venant torsion beam theories are used to derive the governing differential equations of motion and a classical finite element solution is developed.The end moments act about the -axis (lagwise); however, bending in the - plane (lagwise) is not considered and bending occurs in the - plane (flapwise).Thus, the end moments acting in the lagwise direction introduce torsion to the system and create flexural-torsional coupling.Figure 1 illustrates the geometry of the studied system, where , ℎ, and  stand for the beam's length, width, and height, respectively.
The two governing differential equations of the beam are as follows: where  stands for the transverse flexural displacement and  represents the torsional displacement.The derivatives with respect to the length of the beam and time are denoted by a prime (  ) and a dot (⋅), respectively.In (1) and ( 2), the applied moment and force are shown as  and , respectively.The cross-sectional area of the beam is denoted by .The mass density is represented by  and   stands for the polar moment of inertia of the beam.The beam's torsional rigidity () is assumed to be very large compared with its warping rigidity (Γ), and ends are free to warp, that is, state of uniform torsion.As can be observed from (1) and (2), the system is coupled by the end moments, .In order to eliminate the time dependency in (1) and ( 2), simple harmonic vibration is considered and the following transformations are used to describe the transverse and torsional displacements: where  is the circular frequency and  is the time.Ŵ and θ are the transverse and torsional displacement amplitudes, respectively.Upon substituting (3) and ( 4), (1) and ( 2) become The Galerkin method of weighted residuals [40] is employed to develop the integral form of the above equations, where the "hat" signs have been dropped for the sake of simplicity: where  and  (i.e., weighting functions) represent the transverse and torsional virtual displacements, respectively.Performing integration by parts on ( 7) and ( 8) leads to the weak integral form of the governing equations, written as ( ( ( ( ( ( ( ( ( ( (  [ Expressions ( 9) and ( 10) also satisfy the principle of virtual work: where and thus The total virtual work, internal virtual work, and external virtual work are denoted by ,  INT , and  EXT , respectively.The resulting shear force (), bending moment (), and torsional torque (), defined as are zero at the free end and the displacements are set to zero at the fixed boundaries.As a result, the bracketed boundary terms in expressions ( 9) and ( 10) vanish for all boundary conditions.The system is then discretized using elements Nodal DOFs are lateral displacement , rotation (i.e., slope)   , and torsional displacement .The classical finite element formulation is developed using cubic Hermite type polynomial approximations for bending displacement (18) and linear approximations for torsional displacements (19) introduced in the weak integral form of the governing equations such that, for a two-node, three-degree-of-freedom per node element, In ( 18) and ( 19) above  1 and  2 are columns vectors of unknown constant coefficients.The vectors of nodal displacement for bending and torsion are shown below: Thus, where ⟨()⟩ and ⟨()⟩ are both row vectors comprising cubic and linear shape functions for bending and torsion, respectively.The cubic shape functions  1 ,  2 ,  3 , and  4 are The linear shape functions  1 and  2 are defined as This discretizing process leads to the element stiffness and mass and coupling matrices which when assembled together within the FEM code written in MATLAB would result in the linear Eigenvalue problem shown in where  stands for the global stiffness matrix, which is a collection of all the element stiffness matrices.The global mass matrix is symbolized by .Matrix (27a) shown below is the element mass matrix, []  , and matrices (27b) through (27f) are the uncoupled, coupled, and geometric element stiffness matrices.When matrices (27b) through (27f) are assembled together, the final element stiffness matrix would result.This is shown as matrix (27g).Consider where  stands for the element length and  represents the element mass per unit length.The element uncoupled bending stiffness matrix, [  ], is shown below: The final element stiffness matrix is modified due to the presence of the end moment and axial load which contributes the [] geometric matrix, [] torsion matrix, bendingtorsion coupling stiffness matrix, [  ]  , and the torsionbending coupling stiffness matrix, [  ]  .These are added to the bending stiffness matrix, [  ], above, to form the final element stiffness matrix, []  .The geometric and torsion stiffness matrices contributed by the axial load  are shown below: The bending-torsion and torsion-bending coupling stiffness matrices introduced by the end moment  are as follows: Therefore, the final element stiffness matrix, which is a collection of the five submatrices, takes the following form: The solution to the linear Eigenvalue problem in ( 26) is achieved by determining the Eigenvalues and Eigenvectors using a FEM code developed in MATLAB.Various classical boundary conditions are also applied within the MATLAB code.Thus, the natural frequencies and mode shapes of the beam are evaluated.
It is worth noting that if the beam's warping rigidity (Γ) is large compared with its torsional rigidity (), then torsion equation ( 2) changes to a 4th-order differential equation, similar in form to the bending equation (1).While the above FEM formulation is developed for the state of uniform torsion, in that case, following the above-presented procedure and using cubic interpolation functions similar to (24) instead of linear ones (25), to express torsional displacements, the presented formulation can be extended to also include the warping effects.The development of such formulation, however, is beyond the scope of this paper.

Numerical Results
In this section, the validity and practical applicability of the presented FEM procedure are demonstrated through consideration of various examples.A generic beam made of structural steel, subjected to different combinations of axial load, end moment, and end conditions, is first investigated.In a second case study, the comparison is made between the FEM frequency results and limited experimental data available in the open literature.
At first, let us consider a beam made of structural steel ( = 200 GPa and  = 7800 kg/m 3 ), with a length of 8 m, width of 0.4 m, and depth of 0.2 m.The first stage of the numerical tests was to validate the developed FEM code.Due to the lack of analytical results for the preloaded cases, the accuracy of the natural frequency values from the code was established by comparing with the analytical data for an unloaded beam.Table 1 includes the results for the first three natural frequencies for various boundary conditions (with ends free to warp) using the exact method [38] and a 40-element FEM mesh.As it can be observed, the results produced by the exact and the FEM methods are identical and as such the FEM code generates accurate results.
The FEM mesh size was chosen based on a convergence test, as is presented in Figure 3.For this, a cantilever beam, subjected to 1.85 MN of tensile force and 9.21 MN⋅m of end moment, was used.With a focus on higher frequencies, the fifth mode in this case, the FEM convergence was  investigated by increasing the numbers of elements from 10 to 1000, as there were no analytical results available for similar cases.The fifth natural frequency was observed to remain at 88.243 Hz, when 200 or more elements were used.Thus, for the given configuration, the 88.243 Hz value is taken as the exact reference result for the 5th natural frequency.However, since it is apparent that a 40-element mesh results in an error less than 0.01 percent, such a mesh size was considered as reasonable for further studies.Furthermore, it is important to note that convergence of the results for the first four fundamental frequencies was also checked and it was observed that these results converged to the analytical result with an even smaller number of elements.The accuracy of the results produced by the presented FEM method was also verified using a prestressed modal analysis conducted using ANSYS-14 commercial software, where SOLID-187 elements were used.The SOLID-187 element is a higher order, 3D, 10-node element, capable of 6 degrees of freedom (3 translations and 3 rotations) per node.
Tables 2, 3, 4, and 5 present the FEM fundamental frequencies for the beam subjected to different combinations of tensile force and end moment and various boundary conditions, cantilevered (C-F), clamped-clamped (C-C), pinned-pinned (P-P), and pinned-clamped (P-C), respectively.Table 2 also includes a comparison of the results from the developed FEM solution with those from the ANSYS software, which shows a maximum error of 1.3 percent at  = 0 and  = 9.21 MN ⋅ m.These errors may be attributed to the fact that the ANSYS software incorporates the effects  6 and 7, respectively.These tabulated data are also graphically presented in Figure 8. Finally, Figures 9 and 10 depict the bending and torsional components of the  Consider now the axially loaded, industrial aluminum beam ( = 70 GPa and  = 2700 kg/m 3 ), with a length of 1495 mm, and rectangular cross section of 50 mm in width and 10 m min depth, as reported by Laux (2012) [45].The first three natural frequencies of the beam, subjected to different axial loads, are calculated using the presented FEM formulation and are presented in Table 8, along with ANSYS results and experimental data reported in [45].Different mesh sizes are used and, as can be seen from the table, for the first and second natural frequencies, a 5-element course mesh and a 40-element fine mesh both lead to almost identical results.For the 3rd natural frequency, however, there is a slight difference of 0.5% between the results obtained from the two meshes.The FEM and ANSYS results are found to be in excellent agreement with the experimental data, with a maximum difference of less than 1%.For the case in hand, however, there is no experimental buckling data available for comparison (refer to Table 9).

Discussion and Conclusion
The bending-torsion vibration and buckling of beams, subjected to axial load and end moment, were revisited.Neglecting the shear deformation, rotary inertia, and warping effects, the differential equations of motion, coupled by the end moment, were discussed.Exploiting the cubic and linear interpolation functions for bending and torsional displacements, respectively, together with the Galerkin-type weighted residual method, a finite element is developed.Frequency and stability analyses of two illustrative examples are carried out and the FEM results are compared with those obtained from ANSYS software and the data available in the literature.The presented FEM results showed excellent agreement with those obtained from ANSYS and experimental data.As expected, tensile axial load increases the natural frequencies of the beam, indicating an increase in the stiffness of the beam for all classical boundary conditions.When only the end moment is applied, the natural frequencies reduce for all boundary conditions, indicating a reduction in stiffness of the beam.If the end moment is held constant and the tensile load is increased, the natural frequencies increase indicating an increase in the beam stiffness.Conversely, if the tensile load is held constant and the end moment is increased, the beam stiffness reduces.A compressive axial load has the opposite effect and the critical buckling moment reduces with a progressive increase in the compressive load applied.The coupled vibration of the beam, however, is found to be predominantly flexural in the first few natural frequencies (the first three, for the case studied here) and torsion becomes predominant in a higher natural frequency.
Finally, it is worth noting that the presented FEM is designed to account for the axial load, end moment, or combined effects automatically.In contrast, carrying out a similar analysis using the FEM-based software available (ANSYS, e.g.) needs special attention, as the preloads should be introduced in the model through a prestressed analysis.Furthermore, contrary to the analytical formulations, the presented FEM formulation can also be readily extended to the dynamic and stability analyses of (non-)uniform preloaded layered beams and beam structures, to include the warping effects, neglected in the present study, and to analyze beam-like structural components exhibiting other types of coupled behaviour, for example, geometric coupling

Figure 1 :
Figure 1: Beam with axial load and end moment applied at  = 0 and  = .

Figure 2 :
Figure 2: System discretized using elements with 3 degrees of freedom per node.

Figure 3 :
Figure 3: FEM convergence analysis; the fifth natural frequency of preloaded cantilevered beam.

Figure 4 :
Figure 4: Variation of the fundamental frequency with tensile force and end moment; cantilevered (C-F) boundary condition.

Figure 5 :Figure 6 :
Figure 5: Variation of the fundamental frequency with tensile force and end moment; clamped-clamped (C-C) boundary condition.

Figure 7 :
Figure 7: Variation of the fundamental frequency with tensile force and end moment; clamped-pinned (C-P) boundary condition.

Figure 8 :Figure 9 :
Figure 8: Variation of critical buckling compressive force with end moment.

Figure 10 :
Figure 10: Torsional components of the first five natural modes; preloaded cantilevered beam.

Table 1 :
FEM results versus exact data for the first three natural frequencies at  = 0 and  = 0.

Table 2 :
Fundamental frequencies for cantilever beam subjected to force and end moment.

Table 4 :
Fundamental frequencies for pinned-pinned boundary condition (P-P), when force and end moment are applied; a 40element FEM model is used.

Table 6 :
Critical buckling end moment versus compressive force; cantilevered boundary condition.

Table 7 :
Critical buckling compressive force versus end moment; cantilevered boundary condition.