Higher-Order Hierarchical Models for the Free Vibration Analysis of Thin-Walled Beams

This paper addresses a free vibration analysis of thin-walled isotropic beams via higher-order refined theories. The unknown kinematic variables are approximated along the beam cross section as aN-order polynomial expansion, whereN is a free parameter of the formulation. The governing equations are derived via the dynamic version of the Principle of Virtual Displacements and are written in a unified form in terms of a “fundamental nucleus.”This latter does not dependupon order of expansion of the theory over the cross section. Analyses are carried out through a closed form, Navier-type solution. Simply supported, slender, and short beams are investigated. Besides “classical” modes (such as bending and torsion), several highermodes are investigated. Results are assessed toward three-dimensional finite element solutions.Thenumerical investigation shows that the proposedUnified Formulation yields accurate results as long as the appropriate approximation order is considered. The accuracy of the solution depends upon the geometrical parameters of the beam.


Introduction
Many typical aeronautical and space structures involve lightweight, thin-walled, beam-like structures that operate in complex environments.The free vibration characteristics are of fundamental importance in the design of such structures.In particular, the dynamic behaviour of thin-walled beams is richer than that of solid prismatic beams since, besides classical modes (such as bending, torsion, and axial deformation), local higher modes are present.Furthermore, the number of higher modes occurring before the axial one increases as the length-to-side ratio decreases.For these reasons, the accurate modelling and analysis of thin-walled beams (as proposed in this paper) represent an interesting and up-to-date research topic.
A brief overview of recent works about the free vibration of thin-walled beams follows.Matsunaga [1] analysed the natural frequencies and buckling loads of simply supported beams subjected to initial axial forces.Thin rectangular cross sections were investigated.A bidimensional displacement field was assumed.Chen et al. [2] combined the state space method with the differential quadrature method to obtain a semianalytical method for the free vibration analysis of straight isotropic and orthotropic beams with rectangular cross sections.A discussion about properties of the natural frequencies and modes for a Timoshenko beam was presented by van Rensburg and van der Merwe [3].Attarnejad et al. [4] introduced the basic displacement functions that are calculated solving the governing differential equations of transverse motion of Timoshenko's beams by means of the power series method.These functions were applied to the free vibration analysis of nonprismatic beams.Gunda et al. [5] analysed the large amplitude free vibration of Timoshenko beams using a finite element formulation.Transverse shear and rotatory inertia were both accounted for.Different boundary conditions were investigated.Benamar et al. [6] presented a general model for large vibration of thin straight beams.Hamilton's principle was used to obtain a set of nonlinear algebraic equations.Simply supported and clampedclamped boundary conditions were investigated.Tanaka and Bercin [7] studied the natural frequencies of uniform thinwalled beams without cross sectional symmetry.A study on the free vibration of axially loaded slender thin-walled beams was presented by Jun et al. [8].The effects of the warping 2 Mathematical Problems in Engineering stiffness and the axial force were included within Euler-Bernoulli's beam theory.Chen and Hsiao [9] investigated these axial and torsional vibration modes of beams with a Zshaped cross section.The governing equations were derived by the principle of virtual work.The same authors presented in [10] a finite element formulation for the coupled free vibration analysis of thin-walled beams with generic open cross sections.Duan [11] studied the nonlinear free vibrations of thin-walled curved beams with unsymmetric open cross sections via the finite element method.Vörös [12] accounted for the coupling between different vibration modes due to the eccentricity of the cross section as well as to steadystate lateral loads and internal stress resultants.The governing differential equations and boundary conditions were derived using a linearised theory of large rotations and small strains and the principle of virtual work.A finite element model with seven degrees of freedom per node was formulated.Ambrosini [13,14] presented a numerical and experimental study on the natural frequencies of beams with unsymmetric thin-walled open cross sections.Vlasov's theory was modified to include the effects of shear deformation, rotatory inertia, and variable cross section properties.An extension of the previous theory was proposed by de Borbón and Ambrosini [15] to investigate the natural frequencies of thin-walled beams under axial loads.Experimental tests were also carried out in order to verify the proposed theory.Arpaci et al. [16,17] investigated the free vibrations of nonsymmetric thinwalled beams with open cross sections by solving exactly the bending and torsional dynamic equilibrium equations.The coupled bending-torsional behaviour of beams was studied by Banerjee [18] via the dynamic stiffness matrix method.S.-B.Kim and M.-Y.Kim [19] carried out a free vibration as well as a stability analysis of thin-walled tapered beams and frames by means of the finite element method.Zhou-Lian et al. [20] analysed the free vibrations in orthotropic membranes accounting for large deflections.A similar investigation was carried out by by Liu et al. [21] using the L-P perturbation method.
A free vibration analysis of thin-walled isotropic beams is addressed within this paper.The relevance of the proposed analysis is due to the fact that the thinness of the cross section elements results (besides classical bending, torsion, and axial deformation) in local higher modes and corresponding frequencies.The dynamics of thin-walled beams is richer than that of solid prismatic beams and the number of higher modes before the axial one increases as the length-to-side ratio decreases.Therefore, accurate higher-order theories are required.Furthermore, the considered analysis represents a severe test for the proposed models in order to outline both their merits and limitations.In this work, several hierarchical models as well as the classical theories are derived via a Unified Formulation (UF) that has been previously formulated for plates and shells (see Carrera [22] and Carrera and Giunta [23,24]) and lately extended to beams with solid and thin-walled cross sections; see Carrera and Giunta [25], Carrera et al. [26], and Giunta et al. [27].Within this UF, the a priori approximation of the displacement field is written in a compact form.The governing equations are derived through the Principle of Virtual Displacements in terms of a "fundamental nucleus."This nucleus is an invariant of the formulation in the sense that it does not depend upon the theory order of expansion.As a result, a very broad set of beam models that account for transverse shear deformation and cross section in-and out-of-plane warping (although a warping function is not explicitly assumed) can be obtained.Governing differential equations are solved via a Navier-type closed form solution.Open and closed thin-walled simply supported beams are investigated.Slender as well as short simply supported beams are investigated.Open and closed thin-walled cross sections are considered.The modes up to the axial one (which is not included for the sake of brevity) are all considered.The axial mode is disregarded since it is accurately predicted by classical Euler-Bernoulli's and Timoshenko's models.Results are compared with threedimensional finite element models showing that the dynamic behaviour of thin-walled beam structures can be accurately predicted.

Preliminaries
A beam is a structure whose axial extension () is predominant if compared to any other dimension orthogonal to it.The cross section (Ω) is identified by intersecting the beam with planes that are orthogonal to its axis.A Cartesian reference system is adopted: and -axes are two orthogonal directions laying on Ω.The  coordinate is coincident to the axis of the beam.Cross sections that are obtained by the union of  Ω rectangular subdomains with are considered; see Figure 1.Points {(   ,    ) : ,  = 1, 2} are the coordinates of the corner points of a  subdomain.Through the paper, superscript "" represents a cross section subdomain index.The cross section is considered to be constant along .The displacement field is u  (, , ; ) = {  (, , ; )   (, , ; )   (, , ; )} (3) in which   ,   , and   are the displacement components along -, -, and -axis, respectively, and  is the time variable.Superscript "" represents the transposition operator.Stress, , and strain, , vectors are grouped into vectors   ,   that lay on the cross section: and   ,   laying on planes orthogonal to Ω: The strain-displacement geometrical relations are Subscripts "," "," and "," when preceded by comma, represent derivation versus the corresponding spatial coordinate.A compact vectorial notation can be adopted for (6): where D  , D  , and D  are the following differential matrix operators: I being the unit matrix.
Under the hypothesis of linear elastic materials, the generalised Hooke law, holds.The nonzero terms of the material stiffness matrix C are being the Young modulus, ] the Poisson ratio, and  the shear modulus.According to stresses and strain arrangements in ( 4) and ( 5), ( 9) is rewritten as follows: where matrices C  , C  , C  , and C  are

Hierarchical Beam Theories
The variation of the displacement field over the cross section can be postulated a priori.Several displacement-based theories can be formulated on the basis of the following generic kinematic field: u  are unknown functions of the spatial coordinate  and time .  stands for the number of unknowns that depends on the approximation order  (this latter is a free parameter of the formulation).The compact expression in ( 13) is based upon Einstein's notation: an index that is repeated twice implicitly stands for summation.Thanks to this notation, problem's governing differential equations and boundary conditions can be derived in terms of a single "fundamental nucleus."The theoretical formulation is valid for the generic approximation order and approximating functions   (, ).
In this paper, the functions   are the Mac Laurin polynomials.This choice is inspired by the classical beam models.  and   as functions of  can be obtained via Pascal's triangle as shown in Table 1.The actual governing differential equations and boundary conditions due to a fixed approximation order and polynomials type are obtained straightforwardly via summation of the nucleus corresponding to each term of the expansion.According to the previous choice for the polynomial functions, a generic -order displacement field is As far as the first-order approximation is concerned, the kinematic field is Classical models, such as Timoshenko's beam theory (TBT), and Euler-Bernoulli beam theory (EBT), are straightforwardly derived from the first-order approximation model.Within this paper no shear correction coefficient is considered for TBT, since it depends upon several parameters, such as the geometry of the cross section (see Cowper [28] and Murty [29]).Higher-order models yield a more detailed description of the shear mechanics (no shear correction coefficient is required), of the in-and out-ofsection deformations, of the coupling of the spatial directions due to Poisson's effect, and of the torsional mechanics than classical models do.EBT theory neglects them all since it was formulated to describe the bending mechanics.TBT model accounts for constant shear stress and strain components.It should be noted that, for TBT and EBT, the reduced Hooke law for the axial stress/strain relation should be used in accordance with the relaxation procedure classically used to contrast Poisson's locking; see Carrera and Brischetto [30].

Governing Equations
The dynamic version of the Principle of Virtual Displacements (also known as d' Alembert principle) reads where  stands for a virtual variation,   represents the strain energy, and   is the inertial work.

Virtual Variation of the Strain Energy.
According to the grouping of the stress and strain components in (4) and ( 5), the virtual variation of the strain energy is considered as the sum of two contributions: By substitution of the geometrical relations, (7), the material constitutive equations, (11), and the unified hierarchical approximation of the displacements, ( 13), and after integration by parts, (20) reads In a compact vectorial form, The components of the differential stiffness matrix The generic term  gh  (,)  (,) is a cross section moment: that is, a weighted sum (in the continuum) of each elemental cross section area.The weight accounts for material and geometrical distribution.The integral in (24) can be easily computed for the considered domains; see (1) and (2).The following general form for the product   (,)   (,) can be assumed: where   ,   ,   , and   are constant depending upon indexes  and  as in Table 1 and whether differentiation with respect to  and  should be performed or not.The closed form of the generic higher-order moment, therefore, is where the two dots stand for the second derivative versus time and  is the material density.Upon substitution of ( 13), (28) becomes In a compact vectorial form, The components of the matrix M  are where   is Kronecker's delta and For a fixed approximation order, the nucleus has to be expanded versus the indexes  and  in order to obtain the governing equations and the boundary conditions of the desired model.

Closed Form Analytical Solution
The resulting boundary problem is solved via a Navier-type solution.Simply supported beams are, therefore, investigated.
The following harmonic displacement field is adopted: where  is (37) Upon substitution of (35) into (33), the fundamental nucleus of the algebraic eigensystem is obtained: The components of the algebraic stiffness matrix K  are For a fixed approximation order, the eigensystem has to be assembled according to the summation indexes  and .Its solution yields as many eigenvalues and eigenvectors (or modes) as the degrees of freedom of the model.

Numerical Results and Discussion
Beams made of the aluminium alloy 7075-T6 are considered.Mechanical properties are  equal to 71700 MPa and ] equal to 0.3.The material density  is 2.8 ⋅ 10 3 kg/m 3 .Analyses are carried out considering box and I-and C-shaped thin-walled cross sections.The ratio between a representative dimension of the cross section () and the walls' thickness is 20.The cross section sides are such that  = 2 ⋅ 10 −2 m.A lengthto-side ratio equal to 100 and 10 is considered.Slender and relatively short beams are, therefore, investigated.Classical (bending and torsional) and higher modes up to the axial one are considered.The natural frequency of the axial mode is not included in the investigations since it is accurately obtained by classical EBT and TBT.The natural frequencies are put into the following dimensionless form: As far as validation is concerned, results are compared with three-dimensional FEM solutions obtained via the commercial code ANSYS.The quadratic 20-node three-dimensional element "Solid186" is used; see ANSYS v10.0 theory manual [31].The effect of the mesh density is accounted for by considering a fine and a coarse mesh.For the fine mesh, the number of elements along each cross section direction is 40, whereas 20 elements per side are used in the case of the coarse one.The number of elements along the beam axis is such that the element aspect ratio (defined as the ratio between the axial and each cross section side dimension) is equal to ten.The refined three-dimensional solution is addressed as "FEM 3D a " whereas the coarse one is addressed by "FEM 3D b ."The modal shapes of the proposed model and the reference threedimensional solutions are compared via visualisation.Mode visualisation was performed within ANSYS postprocessing environment by imposing at each node the displacement components computed by the proposed models through the ANSYS parametric design language command DNSOL; see Madenci and Guven [32].Since the considered structures are not very complex in terms of geometry, a numerical mode tracking (e.g., by means of the modal assurance criterion, see Allemang [33]) was not necessary.Although the threedimensional FEM solution and the analytical one are different in nature, some considerations about the computational time and effort can be addressed.For the reference FEM simulations, the degrees of freedom (DOFs) over a beam cross section are equal to 2340 (fine solution) and 1140 (coarse solution) in the case of a box cross section.For a fixed approximation order , the DOFs of the proposed solutions are 3( + 1)( + 2)/2 and they do not depend upon the cross section geometry.In the case of the highest considered expansion order ( = 16), they are 459.

Box Cross Section.
A beam with a box cross section as shown in Figure 2 is considered.Several modes were identified before the axial one.Their number increases as the length-to-side parameter decreases.Table 2 presents the natural dimensionless frequencies in the case of / = 100 and  equal to one and two.Since the cross section of the beam has two axes of symmetry, two coincident flexural frequencies are present.Their corresponding modes are indicated as Modes I and II.Torsion is addressed as Mode III.Modes I and II for  = 2 occur before the torsional mode for  = 1.For the flexural modes, the hierarchical beam theories match the FEM solution up to four significant digits with an expansion order as low as four.Classical theories yield good results for the flexural modes where a difference is observed only in the fourth significant digit.The natural frequency related to the torsional mode calls for a higher order of expansion.The error for  = 14 is less than 1%, whereas it is about 4% in the case of a fourthorder model.Consistent with their kinematic hypotheses (the cross section can only translate rigidly on its own plane), classical theories do not yield any torsional mode.The first five dimensionless natural frequencies for the case / = 10 are presented in Table 3. Modes ordering is referred to as the order of apparition in the  = 16 model.The same order is observed in the reference FEM three-dimensional solution, the modes with  being higher than one not considered.Modes I and II are also in this case bending modes.Mode III is a shear deformation on plane  as presented in Figure 3. Mode IV is characterised by bending on an opposite direction of two contiguous sides of the cross section; see Figure 4. Mode V corresponds to torsion.Modes III to V call for high expansion in order.Nevertheless, the shear mode differs from the FEM solution by about 14% even for  = 16.Rather than further increasing the expansion order, the solution accuracy could be improved by using a local approximation approach in a layer-wise sense.This will be matter of future investigations.When the frequencies presented in a table, for a fixed approximation order, are not in an ascending order,  the apparition order is not respected.For instance, Mode V is observed before Mode IV for a ninth-order model.Mode swapping is common in low-order solutions; see Giunta et al. [34], and it is due to a higher cross section stiffness.Classical theories are not suitable when / = 10 because only the two flexural modes are provided.The corresponding values of the frequencies differ from the FEM solution by about four (EBT) and 2% (TBT).

I-Shaped Cross
Section.I-shaped cross section beams are investigated; see Figure 5. Table 4 presents the first three dimensionless natural frequencies for slender beams.The  half-wave number is equal to one and two.The modes are bending on plane , bending on plane , and torsion, respectively.The proposed models all match the FEM solution except for the torsional mode, where the error is about 6% for  = 16.The error for the torsion mode considerably increases when the expansion order decreases.It is about 150%, for instance, for a seventh-order model.In the case of / = 10, six modes are found before the axial one (considering  only equal to one).The first mode is the same as for slender beams.Bending on plane  shifts to the third position, the second mode being a torsional deformation.Modes IV to VI present local deformations ("shell-like" behaviour) of the horizontal and vertical cross section elements.These modes are shown in Figures 6-8.The corresponding frequencies are reported in Table 5.The first three frequencies are accurately predicted.The error for the bending mode on plane  is about five and 8% in the case of TBT and EBT, respectively. = 16 predicts a fourth frequency with a fair accuracy, the error being about 5%.Modes V and VI show the limit of the proposed models and, as a possible solution, a layer-wise local approach towards modelling should be used.

C-Shaped Cross
Section.Finally, beams with a C-shaped cross section are considered; see Figure 9.For the sake of brevity, analyses are carried out only for / equal to 100.
The case of slender beams with  equal to one and two is presented in Table 6.The first three modes, when  is equal to the unit, are bending on plane ,  bending with torsion, and pure torsion.In the case of  = 2, the two bending modes switch the order of apparition.Mode II is pure bending when TBT is used.This results in a sensible error (about 15 and 50% for the two values of , resp.) in the results obtained via classical as well as low-order theories.This example clearly illustrates the importance of higher-order models also in the case of slender beams.The error in the case of the torsional mode for  = 15 is about 2% at worst.

Conclusions
A Unified Formulation has been proposed for the free vibration analysis of thin-walled simply supported beam structures.Via this approach, higher-order models accounting for shear deformations, in-and out-of-plane warping, and rotatory inertia can be straightforwardly formulated.Classical models, such as Euler-Bernoulli's and Timoshenko's, are obtained as particular cases.A closed form, Navier-type solution has been addressed.Beams with thin-walled box as well as I-and C-shaped cross sections have been investigated.Three-dimensional FEM solutions obtained via the commercial code ANSYS have been considered for validation.On the basis of the presented results, it can be concluded that higher-order models are necessary for accurately predicting the frequency of bending modes in C-shaped slender beams.The considered higher-order models yield accurate results in most of the considered cases.One of the leading ideas behind the present paper was to also present the limitation of the proposed models.For this reason, several higher than classical modes have been considered.For some of these cases, few differences from the considered reference threedimensional FEM solutions were observed.This is due to the local "shell-like" deformations of the cross section elements.Besides increasing the cross section expansion order, the adoption of a local layer-wise description of the displacement field can be a possible manner for increasing the solution accuracy.This latter approach will be matter of future work.

) 4 . 2 .
As far as the boundary conditions are concerned, the components of Π  are Virtual Variation of the Inertial Work.The virtual variation of the inertial work is 36) with  ∈ N * representing the half-wave number along the beam axis. = √ −1 is the imaginary unit and  the natural frequency.{  :  = , , } are the maximal amplitudes of the displacement components.The displacement field in (35) satisfies the boundary conditions since  ,  (0) =  ,  () = 0,   (0) =   () = 0,   (0) =   () = 0.
a: 40 × 40 cross section mesh.b: 20 × 20 cross section mesh.c: mode not provided by the theory.
a: 40 × 40 cross section mesh.b: 20 × 40 cross section mesh.c: mode not provided by the theory.