Variable Kinematic One-Dimensional Finite Elements for the Analysis of Rotors Made of Composite Materials

It is well known that composite materials present excellent mechanical properties, such as high specific stiffness and strength, ease of formability, a wide range of operating temperatures, and many others. These properties justify their extensive use in many applications, among which the design of rotors, whose dynamic behavior is worthy of study. For example, interesting experiments were carried out in Refs. [1–3] on either graphite or boron/epoxy shafts. The advantages of orthotropic materials over their isotropic counterparts were pointed out, and useful references were provided for analytical formulations. Despite their limitations, classical beam theories have been extensively used to investigate the critical speeds and instabilities of composite shafts. For instance, with the purpose of proposing an optimization algorithm to define the best lamination scheme, a first shear deformation theory, including the rotatory inertia, was used in Ref. [4]. The optimization was direct in order to maximize the first bending frequency, thus ensuring a sufficient torsional strength by imposing that the lamination angle of a certain number of plies was equal to 45 deg. Later, Bert and Kim developed one-dimensional models, based on Euler–Bernoulli [5] and Bresse–Timoshenko [6] theories, which included bending–twisting coupling. The critical speeds were in good agreement with the results obtained through both shell and experimental approaches. Naturally, when shear effects become important, the Euler–Bernoulli theory does not lead to accurate results. In Ref. [7], Chen and Peng investigated the stability of composite spinning cylinders subjected to compressive loads using the Timoshenko model and the equivalent modulus beam theory (EMBT). The critical speeds and the effects due to a disk location were also studied under these assumptions. A further contribution was made by Chang et al. [8], who developed a simple first-order shear deformation theory for shafts supported by bearings modeled as springs and viscous dampers. The constitutive relations were derived directly from those of three-dimensional continua by means of coordinate transformations.

In order to overcome the limitations of the equivalent modulus approach, Singh and Gupta [9] first proposed a refined beam formulation derived from a layerwise shell theory (LBT) and then, in Ref. [10], presented a modified EMBT that included bending–twisting, shear–normal, and bending–stretching coupling effects. The obtained results were used by Sino et al. [11] to assess their simplified homogenized beam theory (SHBT). In the SHBT, the stiffness parameters are evaluated using an energy formulation that takes into account several contributions, such as Young's modulus, the shear modulus, the distance from the shaft axis, and the thickness of each rotor layer. The results were in good agreement with those presented by Gupta et al. Further attempts have also been made to improve displacement theories. Among these formulations, Librescu et al. provided a higher-order model that incorporates the warping and torsion of thin-walled anisotropic spinning structures. Critical speeds as well as stability were investigated considering thin-walled boxes [12,13] and cylinders [14].

When the assumptions of beam theories are too restrictive, two- and three-dimensional formulations become useful. For instance, in Ref. [15], Kim and Bert compared a number of shell theories (Loo's, Morley's, Love's, Donnell's, and Sanders's theories) to determine the critical speeds of the shaft examined in Ref. [1]. All these models yielded accurate results, except Donnell's theory, which was not effective for long shafts. A further example is Ref. [16], in which Ramezani and Ahmadian combined the layerwise theory and the wave propagation approach for rotating laminated shells under different boundary conditions. Interesting results were presented in Ref. [17], in which the dynamics of rotating cross-ply laminated cylinders was considered. The structure was reinforced with stringers and rings whose stiffness was considered either separately or smeared over the shell surface. As expected, the number of stiffeners had a remarkable effect on the backward and forward frequencies. Moreover, in Ref. [18], a disk–shaft assembly was discretized with shell finite elements, considering a cyclic symmetry of the structure; three analyses were carried out with different numbers of degrees of freedom. This assumption made it possible to limit the computational cost to a great extent. Finally, a comprehensive overview of rotordynamics phenomena with shell approaches can be found in the book by Hua et al. [19].

The present paper has the aim of introducing and comparing a variety of refined one-dimensional beam models in order to study the dynamics of composite rotors. The cross sections of the structures are built with a number of small orthotropic material plates, which are rotated around the rotor–axis of the β angle (see Fig. 1). This approach makes it possible to consider both thin-walled and compact cross sections with symmetric and asymmetric lamination schemes. As far as the thin-walled structures are concerned (cylinders and boxes), higher-order theories are able to detect deformations of the section, which may have an important effect on the overall dynamic behavior. These theories are derived from Carrera's unified formulation (CUF). This procedure makes it possible to consider the order and the type of models as free input parameters. CUF was originally developed for plates and shells [20–22] and was later extended to the beam model [23]. Using one-dimensional finite elements, static [24] as well as free vibration [25] analyses have already been carried out on laminated and sandwich structures, thus demonstrating that refined displacement fields lead to accurate results with a very low computational cost. Recently, Carrera et al. have extended the unified formulation to the study of spinning metallic rotors with thin disks [26] and centrifugally stiffened beams [27].

In light of the encouraging results, the equations of motion (EOM) of spinning composite structures have been written for the first time. The displacement fields are obtained by adopting Taylor-like expansions. The EOM in the rotating frame are solved using the finite element method. In order to assess the new formulation, several linear analyses have been carried out and the related results have been compared with the solutions presented in the literature.

and $r={xP,0,zP}$ is the distance of a generic point P belonging to the cross section from the neutral axis. For sake of clearness, in Appendix  B, the nine components of the fundamental nucleus of the matrix $Kijτs$ are written in an explicit form.

The problem in Eq. (20) is in the classical form and it can be solved with the standard eigensolvers.

This section shows the results related to analyses carried out on composite thin-walled structures. The first case concerned the dynamic study of hollow rotating cylinders whose material and geometrical features are summarized in Table 1. The shaft was discretized with seven four-node beam elements and it was supported by one bearing at each end, whose stiffness k_xx and k_zz was 1740 GN/m. The first critical speeds obtained with the present theories are shown in Table 2 and compared with those found in Ref. [8], adopting the $[90 deg/45 deg/-45 deg/0 deg6/90 deg]$ lamination scheme. The results for both materials are within the interval of reference speeds, and the differences between the three theories (TE2, TE3, and TE4) are essentially negligible. Furthermore, the effects of the aspect ratio for the graphite–epoxy shaft have also been investigated in Table 3. As expected, when the structure becomes thicker and thicker, the critical speeds occur at higher values, and the differences between the three expansions increase; in fact, when the aspect ratio is equal to 2, TE3 and TE4 lead to similar results, whereas TE2 overestimates the critical speed of 10%. In addition, considering the same material, the effects due to the lamination scheme have been evaluated and the results are shown in Table 4. Increasing the lamination angle, with respect to the longitudinal direction, reduces the shaft stiffness and causes a consequent reduction in the critical speeds. In the case in point, some evident discrepancies can be observed between the refined beam models and those found in literature. In fact, when θ is equal to 15 deg, 30 deg, or 45 deg, the present results are closer to speeds obtained with the shell theory SST than those of classical beam models. In order to show the accuracy of the CUF elements, natural frequencies of the boron–epoxy shaft (Table 1), computed with TE elements, are listed in Table 5 and compared with the results found in Ref. [28], where a first-order shear beam theory with a torsion function was conceived, and with those obtained using Di.Qu.M.A.S.P.A.B. [29], a free software based on the differential quadrature method for shells and plates. At least the third-order theory (TE3) is necessary for all the considered cases to obtain accurate results for bending as well as torsional modes. Moreover, it is worth noting that remarkable differences exist between the first natural frequencies of Ref. [28] and those obtained with the other approaches, especially for θ equal to 15 deg, 30 deg, 45 deg, and 60 deg. The use of a constant shear correction factor could be the reason for these important discrepancies, since the shear effects are closely related to the lamination scheme.

In the following numerical illustrations, box beams with thin walls, whose dimensions are shown in Table 1, have been studied. The results are displayed for different cross section height-to-width ratios (⁠$R=b/c$⁠, Fig. 2) and lamination angles (θ), in terms of nondimensional frequency (⁠$f¯=f/fn$⁠), as a function of the speed parameter (⁠$Ω¯=Ω/fn$⁠). The value f_n is the natural frequency of a cantilever graphite–epoxy box beam with R = 1 and θ = 90 obtained with a 2D finite element solution (f_n = 52.649 Hz). Figure 3 shows how the first two frequency ratios vary with the rotational speed for θ equal to 0 deg and 90 deg, using various expansions. As expected, the theory order does not affect the trend of the branches, except that curves start at lower values when the displacement field is enriched, and these differences are more evident for the case in which θ is equal to 0 deg. This fact demonstrates that classical models (and ad-hoc theories!) cannot guarantee the same accuracy for all lamination angles, thus confirming that a simple method to conceive increasingly accurate theories could be a very useful tool. Moreover, in order to evaluate the effect due to the ply angle, the first two frequency ratios are shown in Fig. 4, adopting the TE2 and TE6 expansions. Critical speeds occur at different values, depending on the lamination angle, and maximum and minimum velocities occur for θ = 0 deg and θ = 90 deg, respectively. Similar results were obtained by Librescu et al. in Ref. [12]. In addition, it should be underlined that different lamination schemes determine the occurrence of couplings between different deformation modes. For instance, referring to Fig. 2, two different configurations $[θT/θL/θB/θR]$ were studied: case I$[45 deg/-45 deg/-45 deg/45 deg]$ and case II$[45 deg/-45 deg/45 deg/-45 deg]$⁠. Free vibration analyses were performed to compare the first eight frequencies with those of a 2D finite element solution (Table 6). It should be noted that, for both cases, higher-order theories yield results that are increasingly closer to references. The frequencies (f) and damping (D) ratios were computed as functions of the rotating speed using a number of displacement models (Figs. 5 and 6). Classical theories (EBBT and FSDT, Fig. 5(a)) furnish a qualitatively similar result for the first configuration, but do not predict any instability range. Nevertheless, with the first-order shear deformation theory, a veering of the fourth, fifth, and sixth branches can be detected. Instead, with the TE3 and TE6 expansions, the graphs show dramatic changes, and three instability fields in fact appear within the considered speed interval. This first configuration involves a dominant twist motion and models that are able to describe bending–torsional coupling are therefore needed. In addition, when the expansion order is increased, instability thresholds occur at lower speed values. On the contrary, for the second case, the most important effect is due to shear and, for this reason, the shearable beam model leads to comparable results with those obtained with higher-order theories. For both schemes, the gyroscopic coupling markedly affects the system, causing several veerings of frequency branches.

In the last illustrative example, a rectangular box beam (R = 0.5) was studied with the purpose of considering another kind of anisotropy. Hereafter, the reference frequency f_n is equal to 25.976 Hz. Figure 7(a) has the purpose of showing how the first two frequency ratios vary with the speed parameter when θ is either 0 deg or 15 deg. The analyses were performed adopting four different theories and, for both cases, the curves move toward lower and lower values when the displacement model is enriched. For θ = 0 deg (thin lines), the TE3 branches almost overlap those of the TE6 theory, while discrepancies between these models clearly increase when θ is equal to 15 deg (bold lines). Therefore, in order to ensure a good accuracy for the other lamination schemes, the TE6 expansion was used and the results are shown in Fig. 7(b). It is interesting to note that the curves follow similar paths regardless of the lamination. In fact, after a more or less evident veering of the backward frequency branch, two eigenvalues interact and instability occurs. Furthermore, when θ increases, the instability threshold appears at lower and lower speeds, while the range of the instability field becomes smaller and smaller until it reaches a minimum for θ = 90 deg. The effects of ply angles on the damping are shown in Fig. 7(c).

In the present work, Carrera's unified formulation has been used to study the dynamics of composite rotors. By invoking Hamilton's principle, the equations of motion have been derived and solved with the finite element method. Several kinds of rotors have been considered to assess the new theory and the related results have been compared with those found in literature or with 2D finite element solutions. Higher-order elements were tested on shafts with circular and rectangular cross sections. In light of the results, it is possible to draw the following conclusions:

• The use of finite elements based on higher beam models leads to accurate results regardless to the lamination scheme.

• When the lamination configuration involves bending–torsional coupling, higher theories become mandatory to predict both instability thresholds and critical speeds.

• When composite structures are considered, classical theories yield acceptable qualitative results for a limited range of problems.

Future work could be focused on analyzing other aspects of rotating structures, for instance, by introducing the centrifugal stiffening contribution. The structural models proposed in this paper seem able to describe the dynamics of spinning structures with a high deformability, such as thin-walled cylinders and composite thin disks.