Dynamic Finite Element Analysis of Bending-Torsion Coupled Beams Subjected to Combined Axial Load and End Moment

The dynamic analysis of prestressed, bending-torsion coupled beams is revisited.The axially loaded beam is assumed to be slender, isotropic, homogeneous, and linearly elastic, exhibiting coupled flexural-torsional displacement caused by the end moment. Based on the Euler-Bernoulli bending and St. Venant torsion beam theories, the vibration and stability of such beams are explored. Using the closed-form solutions of the uncoupled portions of the governing equations as the basis functions of approximation space, the dynamic, frequency-dependent, interpolation functions are developed, which are then used in conjunction with the weighted residual method to develop the Dynamic Finite Element (DFE) of the system. Having implemented the DFE in a MATLABbased code, the resulting nonlinear eigenvalue problem is then solved to determine the coupled natural frequencies of illustrative beam examples, subjected to various boundary and load conditions. The proposed method is validated against limited available experimental and analytical data, those obtained from an in-house conventional Finite Element Method (FEM) code and FEMbased commercial software (ANSYS). In comparison with FEM, the DFE exhibits higher convergence rates and in the absence of end moment it produces exact results. Buckling analysis is also carried out to determine the critical end moment and compressive force for various load combinations.


Introduction
Many terrestrial, mechanical, and aerospace structures can be modeled as beams or assemblies of beams, and, therefore, modelling and analysis of such structural elements have been the subject of numerous investigations.Depending on their applications, diverse geometries, loadings, and boundary conditions arise in the structural modeling, leading to a variety of problems.The dynamic, buckling, and vibrational analyses of diverse beam configurations, represented by different geometries and loading scenarios, governed by pertinent theories, have been investigated and reported in the literature.The vibrational analysis of prestressed beams has been the subject of several studies.Neogy and Murthy [1] carried out one of the earliest studies in this area and found first natural frequency of an axially loaded column for two different boundary conditions: pinned-pinned and clampedclamped.Krishn et al. [2], using Rayleigh-Ritz principle, introduced an iterative approximate solution method.Gellert and Gluck [3] investigated the effect of applied axial force on the lateral natural frequencies of a clamped-free beam with transverse restraint.Pilkington and Carr [4] introduced an approximate, noniterative solution for the frequencies of beams subjected to end moment and distributed axial force.Wang et al. [5] used Galerkin's formulation, while Tarnai [6] exploited the more generalized variational technique to investigate the lateral buckling of beams hung at both ends.Later, Jensen and Crawley [7] studied the frequency determination techniques for cases where coupling is caused by warping of composite laminate.They also compared the results of Rayleigh-Ritz and partial Ritz methods with their experimental results.Mohsin and Sadek [8] and Banerjee and Fisher [9] implemented Dynamic Stiffness Matrix (DSM) method to find natural frequencies of an axially loaded coupled beam, while Jun et al. [10] used DSM for vibrational analysis of a composite beam.The DSM method was first introduced by Kalousek [11] for an Euler-Bernoulli beam and ever since has been taken further by many researchers [9,[12][13][14][15][16].
With the advent of more powerful computers in recent years, there has been an increasing interest among researchers to use computational methods in structural stability and vibration analyses.This is mainly due to the fact that the experimental methods are expensive, require extensive testing and measuring techniques, and are limited in their scope of predictions.On the other hand, the analytical solutions are limited to special cases.The classical Finite Element Method (FEM), as the most popular computational technique in solid and structural mechanics, has been extensively utilized by researchers [17][18][19][20].In FEM, fixed shape functions are used to express the field variables in terms of nodal values and to develop the element matrices.Because of their ease of manipulation, Hermite cubic shape functions are commonly used to express elements lateral displacement, resulting in an approximate solution including mass and static stiffness matrices.
In 1998, Hashemi [21] introduced a semianalytical Dynamic Finite Element (DFE) formulation, a hybrid approach that bridged the gap between DSM and FEM methods.Analogous to the conventional FEM, the DFE formulation is based on the general procedure of weighted residual method.However, in contrast to the FEM, the use of frequency-dependent trigonometric shape functions in DFE leads to a frequency-dependent (dynamic) stiffness matrix, which represents both inertia and stiffness properties of the element embedded in a single matrix.As a result, the DFE can be extended to more complex cases, for which a DSM cannot be developed.Since its inception, the DFE method has been extended to vibration analysis of various problems of beamlike structures [22,23].Hashemi and Richard (1999) [22] presented dynamic shape functions and a DFE for the vibration analysis of thin spinning beams.Hashemi and Richard (2000) [23] presented a DFE for the free vibration coupled bendingtorsion beams and investigated the coupled bending-torsion natural frequencies of axially loaded beams and the DFE frequency results were validated against FEM and DSM [9] data and those available in the literature.When compared to the conventional FEM, the DFE generally exhibited much higher convergence rates, especially for higher modes of vibration.
Most of the above-mentioned investigations focused on either uncoupled lateral or geometric/materially coupled stability/vibration of beams.The presence of prestress in structural components can also significantly change the system's stability, dynamic behavior, and response.Helicopter, propeller, compressor and turbine blades, aircraft wing, and rockets internal structure subjected to axial acceleration are some examples of such situations where, at the preliminary design stage, an axially loaded beam model is often used for the dynamic and stability analyses of the system.Also, a beam column with two planes of symmetry, connected through semirigid connections, loaded in the plane of greater bending rigidity by end moments, exhibits coupled torsional-lateral displacements (in the plane of smaller bending rigidity).The former configurations have been thoroughly investigated and reported in open literature.However, studies on the latter cases and the cases of coupled buckling and dynamic behavior of beams subjected to combined axial load and end moment as well as the coupling effects caused by end moments are scarce.Joshi and Suryanarayan [24] developed a closed-form analytical solution for vibrational analysis of a simple uniform beam subjected to both constant end moment and axial load.Later, they unified their solution for different boundary conditions [25] and subsequently developed a general iterative method for coupled flexural-torsional vibration of initially stressed beams [26].More recently, the authors presented a comprehensive study of the coupled bendingtorsion stability and vibration analysis of such elements using the conventional FEM, where the flexural and torsional displacements were expressed using cubic Hermite and linear interpolation functions, respectively [27].
In what follows, a Dynamic Finite Element (DFE) for the coupled flexural-torsional stability and vibration analyses of slender beams, subjected to combined axial force and end moment, is presented.Based on the previously developed applicable governing differential equations of motion, the weighted residual method and integration by parts are exploited to develop the weak integral form of the governing equations.The closed-form solutions of the differential equations governing uncoupled bending and torsional vibrations of an axially loaded beam are used as the basis functions of approximation space to derive the pertinent Dynamic (frequency-dependent) Trigonometric Shape Functions (DTSFs).Introducing the field variables, expressed in terms of the DTSFs and the nodal displacements, into the weak integral form of the governing equation followed by extensive mathematical manipulation leads to the element Dynamic Stiffness Matrix (DSM).The element matrices are then assembled and the boundary conditions are applied to form the system's nonlinear eigenvalue problem, which is finally solved to extract natural frequencies and modes of the system or to evaluate the critical axial load/end moment.It is worth noting that the present DFE is applicable to the members composed of closed sections, with the torsional rigidity, , being very large compared to the warping rigidity, Γ, and with ends free to warp, that is, state of uniform torsion, where the twist rate is constant along the span.However, the presented DFE formulation can also be extended to more complex configurations, such as thinwalled beams with closed or open cross sections, where torsion-related warping effects cannot be neglected.

Theory
Consider a linearly elastic, homogeneous, isotropic slender beam subjected to two equal and opposite end moments, , about -axis and an axial load,  (i.e., loaded in the plane of greater bending rigidity), undergoing linear coupled torsion and lateral vibrations along -axis.Figure 1 depicts the schematic of the problem, where , ℎ, and  stand for the beam's length, width, and height, respectively.Governing differential equations of motion can be developed by defining an infinitesimal element (refer to the authors' earlier work [27]) and by using the following assumptions: (1) The beam is made of linearly elastic material.
(3) The stresses induced are within the limit of proportionality.
(4) The cross section of the beam has two axes of symmetry.
(5) The cross-sectional dimensions of the beam are small compared to the span.
(6) The transverse cross sections of the beam remain plane and normal to the neutral axis during bending, and the beam's torsional rigidity () is assumed to be very large compared with its warping rigidity (Γ), and the ends are free to warp, that is, state of uniform torsion.
The governing differential equations for a prismatic Euler-Bernoulli beam ( and  constant) subjected to static constant axial force () and end moment (  ), undergoing coupled flexural-torsional vibrations caused by end moment, is written as follows [27]: where ()  stands for derivative with respect to  (0 ≤  ≤ ) and ( ̇) denotes derivative with respect to  (time).The internal shear force, (), bending moment, (), and torsional torque, (), are defined as As can be observed from ( 1) and ( 2), the lateral and torsional displacements of the system are coupled by the end moments,   .Exploiting the simple harmonic motion assumption, displacements,  and , are written as where  denotes the frequency and Ŵ and θ are the amplitudes of flexural and torsional displacements, respectively.Substituting ( 4) into ( 1) and ( 2) leads to having Application of the Galerkin weighted residual formulation [19,20] leads to the integral form of the differential equations ( 5) and ( 6), written as where  and  are the weighting functions associated with flexure and torsion, respectively.Performing integrations by parts twice on (7) and once on (8) leads to the following weak integral forms: The above expressions also satisfy the principle of virtual work (PVW), where  is the total virtual work,  INT is the internal virtual work, and  EXT is the external virtual work.For free vibrations,  EXT = 0, and where   and   denote the components of the virtual work associated with flexure and torsion, respectively.It is worth noting that the bracketed boundary terms in ( 6) and ( 7) will vanish regardless of the type of boundary conditions applied, for example, zero displacements and slope,  =   =  = 0, and virtual displacements and slope,  =   =  = 0, at the fixed end (i.e., where the field variables are imposed) and zero resultant internal shear force, (), bending moment, (), and torsional torque, (), at the free end (see expressions ( 3)).
The system is then discretized along the beam span by a certain number of 2-noded, 6-DOF elements, such that (see Figure 2) where the discretized (element) weak form equations are with the geometric and material parameters all assumed to be constant per element.At this point, a conventional FEM can be developed, where polynomial interpolation functions are used to express the field variables in terms of nodal variables (see, e.g., [19,27]).The Dynamic Finite Element (DFE), as mentioned earlier in introduction section, is a highly convergent and efficient hybrid approach combining the classic Finite Element (FEM) and Dynamic Stiffness Matrix (DSM) methods.The DFE method is equipped with frequency-dependent trigonometric shape functions inspired by DSM formulation with averaged value parameters over each element.The solutions of uncoupled governing differential equations are used as basis functions of approximation space to derive the dynamic shape functions, which are then exploited to find the element (frequency-dependent) dynamic stiffness matrix.
The DFE formulation starts with the weak form of element equations ( 13), where further integration by parts is applied on the first two terms of each equation, leading to the following form of the equations: Substituting  = / (0 ≤  ≤ 1) in the above equations results in the following element nondimensional equations: The closed-form solutions of the integral terms ( * ) and ( * * ), respectively, can be written as [9]  () =  1 sin () +  2 cos () +  3 sinh () where  1,2,3,4 and  1,2 are constants and with Specific linear combinations of components of solutions ( 16) are used as the basis functions of approximation space, and the general FEM procedure is exploited to derive the interpolation functions used to express element variables in terms of the nodal properties (which also satisfy the integral terms marked as ( * ) and ( * * )).For this purpose, the nonnodal solution functions,  and , and the test functions,  and , are written in terms of generalized parameters ⟨⟩, ⟨⟩, ⟨⟩, and ⟨⟩ as follows: with the dynamic (frequency-dependent) basis functions of approximation space defined as Replacing the generalized parameters, ⟨⟩, ⟨⟩, ⟨⟩, and ⟨⟩ with the nodal variables, and ⟨ 1  2 ⟩, respectively, expressions (20) can be rewritten as where the matrices [  ]  and [  ]  are defined as By combining expressions (20) and above matrices, [  ]  (23), and [  ]  (24), the nodal approximations for flexural displacement, (), and torsional displacement, (), are written as where ⟨()⟩  and ⟨()⟩  are the frequency-dependent (dynamic) trigonometric shape functions for flexure and torsion, respectively.Then, (25) can also be rewritten as Shock and Vibration where, The expressions of flexural frequency-dependent trigonometric shape functions are as follows: where and the torsional trigonometric shape functions are written as where   = sin() [21].
Using expressions (15) and the shape functions ( 28 ] ] ] and the two coupled element matrices are written as The element dynamic stiffness matrix, [()]  , is determined by adding these six coupled and uncoupled submatrices and the global dynamic stiffness matrix, [()], is then obtained by assembling all the element matrices and applying the system boundary conditions, carried out using a MATLAB code.According to the principle of virtual work (10), for arbitrary virtual displacement, {  }, the nonlinear eigenvalue problem of the system resulting from the DFE process is written as follows: where {  } contains all the nodal displacements of the system.The natural frequencies of the system would be the nontrivial solutions of expression (37), that is, values of  which yields a zero determinant for the global dynamic stiffness matrix, |()| = 0. Consequently, the system natural modes, {  }, can be found by extracting data from corresponding eigenvectors.
It is worth noting that the basis functions of approximation space have been specifically designed such that when the frequency of oscillations tends to zero, the roots, , , and  of the characteristic equations also tend to zero.Consequently, the limits of basis function (17) and (18), formed as linear combinations of the closed-form solutions to the starred characteristic equations in expressions (15), change to those used in the classical beam FEM.In other words, expression ( 17) and ( 18) become cubic Hermite and linear polynomials, respectively, commonly used as flexural and torsional basis functions in FEM.Subsequently, element dynamic stiffness matrices (31) through (36) reduce to the static stiffness matrices obtained from conventional FEM; that is, DFE transformed to FEM (see, e.g., [27]).However, this would not be possible if the individual components of solution functions (21) were simply taken as the basis functions.
It is also 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 (8).While the presented DFE is developed for the state of uniform torsion, the formulation can be extended to include the warping effects.This could be done by following the abovepresented procedure and using the four dynamic interpolation functions in a form similar to ( 28)-(29) instead of (30) currently used to express torsional displacements.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 DFE method are demonstrated through various illustrative examples.At first, the DFE was validated against limited existing analytical data in the literature.A generic beam made of structural steel, subjected to different combinations of axial load, end moment, and various end conditions, was then investigated.The DFE results were compared with those obtained from an in-house FEM-based code [27] as well as ANSYS and limited experimental data available in the open literature.
Consider a slender beam of rectangular cross-sectional area, with Young modulus  = 200 GPa, density  = 7800 kg/m 3 (steel), length of 8 m, width of 0.4 m, and depth of 0.2 m.The DFE convergence for the beam's 5th natural frequency was first established, as depicted in Figure 3.The relative error is found based on the reference values found from analytical closed-form solution ( = 1.23 MN,   = 0) presented by Joshi and Suryanarayan [24].A comparison between the DFE method and conventional FEM with regard to the convergence rate is illustrated in Figure 4.For the 5th natural frequency, a 5-element DFE model produces results with less than 0.2 percent error compared to the analytical result or an error less than 0.1 percent using an 8element mesh, whereas an 8-element FEM model shows a 0.3 percent error, that is, three times larger error than DFE.The convergence tests were also carried out for the 1st through 4th frequencies and showed even better performances for DFE over conventional FEM.Tables 1-4 present the DFE and FEM results for the beam's fundamental frequency, subjected to various combinations of preloads, and different boundary conditions: cantilever (C-F), clamped-clamped (C-C), pinned-pinned (P-P), and pinned-clampled (P-C).Preload combinations include two end moments of   = 0 and 6.14 (MN⋅m) and four axial loads of  = 0, 0.62, 1.23, and 1.85 MN.
It is worth noting that, in absence of end moment, the equations of motion are uncoupled, and, as a result, the DFE method yields exact results [21].In other words, when   = 0, the DFE method leads to exact frequency data.However, when the system is subjected to an end moment (  ̸ = 0), an exact solution does not exist and, therefore, the frequencies  1, the relative errors for the beam's fundamental frequency (subjected to  = 0.62 MN and   = 6.14 MN⋅m) are found to be 0.11% and 3.67% for 5-element DFE and FEM models, respectively.In order to investigate the effect of axial force and end moment on the beam stability, a buckling analysis is carried out.Table 5 presents the beam's critical buckling moments for different applied axial forces, and Table 6 shows the buckling forces for different applied end moments (  ), obtained using a 5-element DFE model.
Figures 5-8 are graphical representations of the results presented in Tables 1 through 4.These figures illustrate the variation of the system's fundamental frequency with tensile axial force () and end moment.
Figure 9 illustrates variation of critical buckling end moment ( -cr ) with axial force () for cantilevered boundary condition, and Figure 10 depicts critical buckling compressive force ( cr ) with end moment (  ) for cantilevered boundary condition.Figures 11 and 12 show bending and torsion components of the mode shapes, respectively, for a cantilevered system subjected to a tensile force of  = 1.85 MN and end moment of   = 9.21 MN⋅m.
For further validation, the DFE and conventional FEM (in-house code and ANSYS) methods were used to reproduce the first three flexural natural frequencies of an axially loaded clamped-clamped beam, for which the experimentally evaluated values are available in the open literature [28].It is worth noting that, in this case, the end moment is null,   = 0, and as a result the system exhibits uncoupled displacements.Therefore, as discussed before, the DFE frequency results are  exact, even if only one element is used (see also [21,22]).In this experiment, the beam is made of Aluminum with  = 2700 kg/m 3 ,  = 26 GPa,  = 70 GPa, and dimension of beam is  ×  ×  = 1290 × 75 × 35.As can be seen from the results shown in Table 7, although there is a very close match between the frequencies obtained from all the methods, the DFE results match the best experimental data, in general, and for higher modes in particular.Additionally, as it can be observed from Table 7, the frequency results obtained from in-house and commercial FEM codes are   virtually identical, as they both use the same number of elements and interpolation functions.

Discussion and Concluding Remarks
A Dynamic Finite Element (DFE) formulation was developed to conduct vibrational and buckling analysis of beams subjected to axial load and end moment.Neglecting the shear deformation, rotary inertia, and warping effects, the bendingtorsion coupled nature of vibration caused by applied end moment about the perpendicular axis to bending axis was demonstrated.The DFE results were compared with exact results from DSM, those obtained from FEM and ANSYS models, and limited existing experimental results for the simple case of an axially loaded beam.For the last case, contrary to FEM, the results obtained from a one-element DFE model are exact within the limits of the theory.Similar comparisons were made between DFE, FEM, and ANSYS results for an axially loaded beam subjected to various end moments, and higher rates of convergence for DFE were observed for various classical boundary conditions.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 subjected to end moment only, the natural frequencies decrease for all boundary conditions, indicating a reduction in stiffness of the beam.Holding the end moment constant, the natural frequencies directly increase with tensile load (increased beam stiffness).A compressive axial load has the opposite effect and the critical buckling moment reduces with a progressive increase in the compressive load applied.Conversely, if the tensile load is held constant, the beam stiffness is inversely related to the magnitude of end moment.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 at a higher natural frequency.In summary, the DFE is shown to be a powerful method, exhibiting higher convergence rates than FEM in free vibrational and stability analyses of preloaded beams, in general, and particularly when the higher modes of vibration are of interest.It is also worth noting that the presented DFE is designed to account for the axial load, end moment, or combined effects automatically.Carrying out a similar analysis using the FEM-based software available (e.g., ANSYS) needs special care, as the preloads should be introduced in the model through a prestressed analysis.Furthermore, unlike the analytical formulations, the presented DFE can also be readily extended to the dynamic and stability analyses of more complex cases, such as nonuniform preloaded layered beams and beam structures, and include the warping effects neglected in the present study.

Figure 1 :
Figure 1: Schematic and coordinate system of the problem, with axial load and end moment applied at  = 0 and  = .

Figure 2 :
Figure 2: Discretized domain along the beam span.

Figure 5 :
Figure 5: Variation of natural frequencies when tensile force and end moment are applied for cantilevered boundary condition.

Figure 6 :
Figure 6: Variation of natural frequencies when tensile force and end moment are applied for clamped-clamped boundary condition.

Figure 7 :
Figure 7: Variation of natural frequencies when tensile force and end moment are applied for pinned-pinned boundary condition.

Figure 8 : 5
Figure 8: Variation of natural frequencies when tensile force and end moment are applied for pinned-clamped boundary condition.

Figure 9 :
Figure 9: Variation of critical buckling end moment ( -cr ) with axial force () for cantilevered boundary condition.

Table 5 :
Critical buckling moment for cantilevered boundary condition with varying compressive force (5-element DFE model).

Table 7 :
The first three bending natural frequencies (Hz) of the clamped-clamped aluminium beam for various axial loads: 1-element DFE, 5-element FEM, and ANSYS models.