Experimental analysis of the SHBT approach for the dynamic modeling of a composite hollow shaft

Journal Pre-proofs Experimental analysis of the SHBT approach for the dynamic modeling of a composite hollow shaft Paulo C.P.F. Barbosa, Vergílio T.S. Del Claro, Marcelo S. Sousa Jr., Aldemir A. Cavalini Jr., Valder Steffen Jr. PII: DOI: Reference:

S0263-8223(18)34426-X https://doi.org/10.1016/j.compstruct.2020.111892 COST 111892

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

6 December 2018 10 October 2019 7 January 2020

Please cite this article as: Barbosa, P.C.P., Del Claro, V.T.S., Sousa, M.S. Jr., Cavalini, A.A. Jr., Steffen, V. Jr., Experimental analysis of the SHBT approach for the dynamic modeling of a composite hollow shaft, Composite Structures (2020), doi: https://doi.org/10.1016/j.compstruct.2020.111892

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier Ltd.

Experimental analysis of the SHBT approach for the dynamic modeling of a composite hollow shaft Paulo C. P. F. Barbosa, Vergílio T. S. Del Claro, Marcelo S. Sousa Jr., Aldemir A. Cavalini Jr.*, Valder Steffen Jr. LMEst – Structural Mechanics Laboratory, Federal University of Uberlândia, School of Mechanical Engineering, João Naves de Ávila Avenue, 2121, Building 1O, Uberlândia, MG, 38408-144, Brazil * Corresponding author: [email protected] Abstract Composite materials have been extensively used in engineering applications over the past few decades, enabling new possibilities for engineering design. The limitations faced by composites have been eliminated through the development of new resins, fibers, and configurations. Thus, composite shafts appear to be an interesting solution for many rotating machinery applications due to their associated low weight, high-strength, high-temperature resistance, mechanical properties manipulation, and thermal cycling resistance as compared with metallic shafts. For the modeling of composite shafts, simplified theories are often used. In the present contribution, a finite element model (FE model) based on the Simplified Homogenized Beam Theory (SHBT) was used to predict the dynamic behavior of a thick-walled composite hollow shaft. The FE model of the shaft is based on the Timoshenko beam theory. A test rig presenting a thick-walled carbon-epoxy composite hollow shaft, two aluminum discs, and two self-aligning ball bearings was used to validate the mentioned FE model. In this case, numerical and experimental frequency response functions (FRFs) and unbalance responses were compared. Results show promising outcomes from the considered FE model regarding the prediction of the first critical speed and the corresponding unbalance responses of the composite rotor. Keywords: composite hollow shaft, rotating machines, SHBT approach, finite element model, dynamic behavior. 1 Introduction The use of composite materials has been steadily growing over the years, mostly in naval, aeronautical, and automotive industries, due to the range of possibilities in obtaining suitable characteristics for specific applications. Rotor dynamics is one of the areas with a great interest in the use of composite shafts since they are a viable solution to overcome the intrinsic limitations of metallic shafts [20]. For systems operating under subcritical conditions (rigid rotors), the goal is to decrease weight and maximize the torque transmission. The low weight of composite shafts allows for faster acceleration and deceleration as compared with conventional rotating machines [4]. It also greatly reduces effects related to rotational inertia, which are undesired in many applications. In supercritical operations

(flexible rotors), the vibration responses associated with the shaft bending, dynamic stress, stability, and fatigue should be carefully evaluated [13]. A forthcoming application for composite shafts has to do with on onboard systems due to the inherent low weight of composite materials. In onboard conditions, the bearings supporting the rotor is subjected to external base excitations [25]. Onboard conditions for rotating machines are observed on high-sea or offshore platforms, vehicles in general, earthquake-resistant machinery, etc. These systems can greatly benefit from the application of composite shafts, regarding all the previously presented advantages, being the most important one its low weight, which reduces the inertia, the undesirable dynamic effects and subsequent instabilities [9]. In composite material shafts, it is possible to change stiffness and damping properties by manipulating some characteristics, such as fiber and matrix composition, fiber orientation, number of layers, stacking sequence, layer thickness, and other geometric properties [3]. This allows for critical velocities to be conveniently changed according to the project requirements for the rotor system. Additionally, it is possible to attenuate the vibration amplitudes when the system undergoes critical speeds [20], depending on the structural damping properties of the composite [16]. The internal structural damping can significantly affect the dynamic behavior of composite shafts. The composite material can be used to reduce the vibration amplitude of the shaft at critical speeds, however, it may also make the system unstable [5, 20]. In rotors with metallic shafts, the internal damping influence can be omitted in most cases. Nevertheless, in composite shafts, this influence can be up to twice as large as on a conventional shaft [29]. In this context, the characterization of the internal damping is important to the design of rotating machines with composite shafts aiming to establish safe operating conditions [8, 23]. Simplifying hypotheses are commonly used to model composite shafts, allowing for the dynamic behavior of the system to be sufficiently accurate represented, whilst maintaining a reasonably low computational cost. Various recent FE formulations based on the homogeneous beam theory, equivalent layer theories, and the costlier layerwise models have been proposed for the analysis of composite materials [1, 2, 17, 18, 30]. Nevertheless, the logical reasoning behind these methodologies is to produce a model complex enough to properly represent the desired system, with a computational cost that is low enough for the problem to be properly solved. The Equivalent Modulus Beam Theory (EMBT) [21, 27] can be effectively used for the modeling of hollow composite shafts and even extended to represent rotating dynamics of machines. However, this EMBT model has some limitations, mainly regarding the coupling of the vertical and horizontal directions found in asymmetric composite shafts [12]. The Layerwise Beam Theory (LBT) [23] presents a satisfactory prediction of the composite shaft dynamics in rotating system applications, but have an exponentially higher computational cost, given the complexity of the beam element utilized. LBT is based on the Layerwise Shell Theory (LST), considering that the cross-section of the beam suffers no distortion whatsoever [21]. Even alternative solutions (i.e., Rayleigh-Ritz method) were recently used in the modeling of composite shafts [8]. Intermediate approaches were evaluated by researchers, resulting in the EMBT model [3] and the Simplified Homogenized Beam Theory (SHBT) [24]. In this study, the SHBT model [24] was used to represent the dynamic behavior and unbalance responses of an experimental test rig with a thick-walled carbon fiber composite hollow shaft. The rotor system was evaluated as operating under sub-critical rotation speeds. As the SHBT model is originally intended for thin-walled rather than

thick-walled shafts, its accuracy depends on the relation between wall thickness and external radius. Regardless of the SHBT model apparent simplicity, it accounts for the internal damping of the composite. The SHBT method associated with the Timoshenko beam theory for an FE model representation allows for modeling any stacking sequence of the plies, considers the structural damping of the composite, and also takes into account the shear stress acting transversely to the shaft. Thus, the main objective of this contribution is to evaluate numerically and experimentally an FE model that is more accurate and precise than other alternative approaches. Aiming at presenting novel experimental results and providing a proper verification of the considered methodology, an experimental validation was performed. Thus, an experimental test rig composed of a thick-walled carbon-epoxy composite hollow shaft, two aluminum discs, and two self-aligning ball bearings was used. Various experimental measurements were carried out, acquiring displacement data both along the vertical and horizontal directions of the two discs. In order to be sufficiently complex for the validation of the FE model, a thick-walled carbon-epoxy shaft with 20 layers (wall thickness to radius ratio of 32.1%) was used. The numerical results are presented both in time and frequency domains, overlaid with experimental results as appropriate, demonstrating the ability of the SHBT model to represent the dynamic behavior of composite shaft rotors. Some physical parameters that depend on the curing process and are especially difficult to be experimentally determined with accuracy, such as the Young’s and transverse stress modules for each individual ply of the composite, they were adjusted via an optimization procedure [15] to ensure a good agreement between the model and experimental results. For this aim, the model updating procedure was based on numerical and experimental FRFs. 2 Rotor model The formulation associated with the rotor FE model and the physical considerations made for its development are hereby detailed. This section will be divided accordingly, aiming at describing the considered SHBT FE model. 2.1 Equations of motion Equation (1) shows the differential equation that represents the dynamic behavior of a flexible metallic rotor operating on a steady state condition [14, 3]. 𝑴𝛿 + [𝑫 + Ω𝑫𝑔]𝛿 + 𝑲𝛿 = 𝑾 + 𝑭𝑢

(1)

where 𝑴 is the mass matrix, 𝑫 is the damping matrix associated with the bearings, 𝑫𝑔 represents the gyroscopic effect, and 𝑲 is the stiffness matrix. The vector δ represents the generalized displacements of the shaft (lateral displacements) and Ω is the rotational speed. 𝑾 stands for the weight of the rotating parts and 𝑭𝑢 represents the unbalance forces. Considering the dissipative effects associated with composite materials [24], Eq. (1) is modified as given by Eq. (2). 𝑴𝛿 + [𝑫 + Ω𝑫𝑔 + 𝑫𝑖]𝛿 + [𝑲 + Ω𝑲𝑖]𝛿 = 𝑾 + 𝑭𝑢

(2)

so that 𝑫𝑖 and 𝑲𝑖 are the internal damping and the stiffness matrices, respectively, both associated with the composite material. This formulation for the composite considers a FE model based on the Timoshenko beam model [3, 19]. In this formulation, the shaft is characterized by its strain and kinetic energies, the discs and unbalance mass are determined from their kinetic energy only, and the bearings are taken into account through the virtual work produced by their supporting forces. Therefore, to obtain the equations of motion for the rotor (Eq. (2)) one must proceed to the Lagrange formalism. The kinetic and potential energies for the composite shaft are presented next. The formulation of the kinetic energies associated with the discs and unbalance mass and the virtual work produced by the bearings supporting forces are described in [14]. 2.2 Composite shaft kinetic and potential energies The Lagrange model for the composite shaft rotor is given by:

( )―

𝑑 ∂𝑇 𝑑𝑡 ∂𝒒𝑖

∂𝑇 ∂𝒒𝑖

∂𝑈

(3)

+ ∂𝒒𝑖 = 𝑄𝑖

where 𝑇 and 𝑈 are the kinetic and potential energies, respectively, and 𝑄𝑖 is the vector of generalized efforts. Concerning the shaft contribution to the overall dynamics, a modified beam element based on the classical Timoshenko beam element is adopted. This element is able to represent the coupling between the directions 𝑢 and 𝑤 (see Fig. 1), presenting 4 degrees of freedom (DoFs) per node, being 2 linear displacements and 2 rotations, resulting 8 DoFs per element [7].

Figure 1: Composite shaft model for a single element. Formulating the shape functions for this element [14], the rotation 𝜃 is the partial derivative of 𝑤 relative to 𝑦, and in a similar way, the rotation 𝜑 is the negative partial derivative of 𝑢 relative to 𝑦. The transverse strain fields for each element are described by a third-degree polynomial, as given by Eqs. (4) and (5).

{𝑤𝑢 == 𝑵𝑵 𝒒𝒒 , 𝑤ℎ𝑒𝑟𝑒 {𝒒𝒒 == {{𝑢𝑤 ,𝜑,𝜃 ,𝑢,𝑤 ,𝜑,𝜃 }} 1 𝑢

𝑢

1

1

2 𝑤

𝑤

1

1

2

2

2

2

(4)

{

{ = {1 ―

𝑵1 = 1 ― 𝑵2

3𝑦2 2

𝐿 3𝑦2 2

𝐿

+ +

2𝑦3 3

, ―𝑦+

𝐿 2𝑦3 3

𝐿

2

, 𝑦―

2𝑦 𝐿

2𝑦2 𝐿

𝑦3

3𝑦2

― 𝐿2, 3

𝑦

― 𝐿2,

2

𝐿 2

3𝑦

2

𝐿

―

2𝑦3 3

𝐿

,

3

―

2𝑦

3

𝐿

, ―

} }

𝑦2 𝐿

― 𝐿2

2

3

𝑦 𝐿

𝑦3

𝑦

+ 𝐿2

(5)

being 𝑵1 and 𝑵2 the shape functions of the beam element under bending, and L is the element length (0 ≤ y ≤ L). The shaft kinetic energy Ts can be determined by Eq. (6), as follows: 𝜌𝑆 𝐿

𝜌𝐼 𝐿

𝐿

𝑇𝑆 = 2 ∫0(𝑢2 + 𝑤2)𝑑𝑦 + 2 ∫0(𝜑2 + 𝜃2)𝑑𝑦 +2𝜌𝐼Ω2 +2𝜌𝐼Ω∫0𝜑𝜃𝑑𝑦

(6)

being 𝜌 the density of the composite shaft, 𝑆 the transverse section area, and 𝐼 the inertia of the transverse section area. By substituting Eqs. (4) and (5) into Eq. (6), one obtains the compact form of the shaft kinetic energy. Subsequently, by applying this compact form to the Lagrange equation, Eq. (3), one obtains Eq. (7). 𝑑 ∂𝑇𝑆 𝑑𝑡 ∂𝒒

( )―

∂𝑇𝑆 ∂𝒒

= (𝑴𝑆 + 𝑴𝑇)𝒒 +Ω𝑫𝑠𝒒

(7)

where the boldfaced terms represent the 8x8 FE matrices for the shaft. More precisely, 𝑴𝑆 is the shaft mass matrix, 𝑴𝑇 includes the shaft rotation inertia secondary effects, and 𝑫𝑠 is the shaft gyroscopic effect matrix. Later, matrix 𝑴 will represent the sum 𝑴𝑆 + 𝑴𝑇. The shaft strain energy is obtained from the analysis of the terms of strain 𝜀 and stress 𝜎. Assuming that the shaft itself is symmetric and 𝐶 is the shaft geometric center, one can characterize the physical model shown in Fig. 2.

Figure 2: Physical model. Adapted from [14]. In this case, 𝑢 and 𝑤 are the displacements in inertial coordinates with respect to the 𝑋 and 𝑍 directions, respectively; 𝑢 ∗ and 𝑤 ∗ are their equivalent rotating coordinate with respect to the axes 𝑥 and 𝑧. One can therefore determine the longitudinal strain of an arbitrary point 𝐵, as given by Eq. (8).

∂2𝑢 ∗

∂2𝑤 ∗

𝜀 = ― 𝑥 ∂𝑦2 ― 𝑧 ∂𝑦2 ,

𝑤ℎ𝑒𝑟𝑒:

{𝑤𝑢

∗ ∗

= 𝑢cos (Ω𝑡) ― 𝑤sin (𝛺𝑡) = 𝑢sin (Ω𝑡) + 𝑤cos (𝛺𝑡)

(8)

This procedure leads to the shaft strain energy U, disregarding axial stresses (either external or induced by coupling with the transversal deformation of the shaft), as shown by Eq. (9). 1

𝑈 = 2∫𝑉𝜀𝑡𝜎𝑑𝑉.

(9)

Considering Hook’s law (𝜎 = 𝐸𝜀), one can modify Eq. (9), substitute Eq. (8) into the modified equation, and expand it, thus resulting Eq. (10).

𝐸 𝐿

[

∂2𝑢 ∗ 2

( )

𝑈 = 2∫0∫𝑆 𝑥2

∂𝑦2

∂2𝑤 ∗ 2

( )

+ 𝑧2

∂𝑦2

∂2𝑢 ∗ ∂2𝑤 ∗

+ 2𝑥𝑧 ∂𝑦2

∂𝑦2

]𝑑𝑆𝑑𝑦

(10)

Then, by considering that the shaft cross-section remains circular, one can reduce Eq. (10) to Eq. (11). 𝐸𝐼 𝐿

[

𝑑2𝑵𝑡1𝑑2𝑵𝑡1

𝑈 = 2 ∫0 𝒒𝑡𝑢 𝑑𝑦2

2

𝑑𝑦

𝑑2𝑵𝑡2𝑑2𝑵𝑡2

𝒒𝑢 + 𝒒𝑡𝑤 𝑑𝑦2

𝑑𝑦2

]

𝒒𝑤 𝑑𝑦

(11)

where 𝐸 is the shaft Young’s modulus, 𝒒𝑡𝑢 represents the vector of DoFs along the plane XY (see Fig. 1), and 𝒒𝑡𝑤 represents the vector of DoFs along the plane YZ. By applying the Lagrange equation (Eq. (3)) into Eq. (11), one can obtain the stiffness matrix 𝑲, as presented by Eq. (12). This equation considers the shear stresses corresponding to the Timoshenko beam theory. ∂𝑈 ∂𝒒

= 𝑲𝒒

(12)

where 𝒒 contains all DoFs of the beam element, as illustrated in Fig. 1. 2.3 Internal damping The internal damping found in composite shafts is treated as an equivalent viscous damping. It is mathematically represented by using the so-called Kelvin-Voigt rheological model [24], as shown in Fig. 3.

Figure 3: Kelvin-Voight element representation. Adapted from [28].

The representative element of this approach is a spring-damper system, associated in parallel, so that the spring is the elastic fraction of the deformation, following Hooke’s law, and the damper is assumed to be linear, with a resulting stress expressed as a function of the strain rate. The relationship between the stress 𝜎 and the deformation 𝜀 of the Kelvin-Voigt [28] model is given by Eq. (13).

{σε == σε ++ εσ

σ = 𝐸𝜀 + 𝜂𝜀, 𝑤ℎ𝑒𝑟𝑒:

1

2

1

2

(13)

where 𝜂 is the viscosity and 𝜀 is the deformation rate. The relaxation time 𝛼 is given by 𝛼 = 𝜂 𝐸. This is a characteristic associated with the composite material structural damping. Thus, the internal dissipative effects of the composite are described by the virtual work 𝛿𝑊 as given by Eq. (14).

[

𝐿

∂2𝑢 ∗ ∂2𝛿𝑢 ∗

𝛿𝑊 = ∫0∫𝑆𝛼𝐸 𝑥2 ∂𝑦2

2

∂𝑦

∂2𝑢 ∗ ∂2𝛿𝑤 ∗

+ 𝑥𝑧 ∂𝑦2

2

∂𝑦

∂2𝑤 ∗ ∂2𝛿𝑢 ∗

+ 𝑥𝑧 ∂𝑦2

2

∂𝑦

]𝑑𝑆𝑑𝑦

(14)

{𝑤𝑢 == 𝑵𝑵 𝒒𝒒 .

(15)

∂2𝑤 ∗ ∂2𝛿𝑤 ∗

+ 𝑥2 ∂𝑦2

∂𝑦2

in which

{𝑤𝑢

∗ ∗

= 𝑢cos (Ωt) ― 𝑢Ωsin (Ωt) ― 𝑤sin (Ωt) ― 𝑤Ωcos (Ωt) , 𝑤ℎ𝑒𝑟𝑒: = 𝑢sin (Ωt) + 𝑢Ωcos (Ωt) + 𝑤cos (Ωt) + 𝑤Ωsin (Ωt)

1 𝑢

2 𝑤

Applying Eq. (15) into Eq. (14), one can obtain Eq. (16).

𝐿

[(

𝐹𝑖 = ―𝛼𝐸𝐼∫0

𝛿𝑊 = ― 𝐹𝑡𝑖𝛿𝒒, 𝑤ℎ𝑒𝑟𝑒: 2

𝑑

𝑵𝑡1𝑑2𝑵1 2

2

𝑑𝑦 𝑑𝑦

+

𝑑2𝑵𝑡2𝑑2𝑵2 2

2

𝑑𝑦 𝑑𝑦

)𝒒 + Ω (

𝑑2𝑵𝑡2𝑑2𝑵1 2

2

𝑑𝑦 𝑑𝑦

―

)𝒒]𝑑𝑦

𝑑2𝑵𝑡1𝑑2𝑵2 𝑑𝑦2 𝑑𝑦2

(16)

being 𝑫𝑖 and 𝑲𝑖 of Eq. (2) obtained from the manipulation of Eq. (16) to determine its matrix form. 2.4 Homogenization method In the present contribution, the SHBT theory was used to model the composite shaft. The homogenization method considers equivalent properties for each beam section. Nevertheless, this method demands a first coordinate system fixed to the shaft and following the inertial coordinates, and a second coordinate system also coupled to the shaft, but following the fibers orientation for each orthotropic composite ply, as can be seen in Fig. 4. For the present model, a general transverse orthotropic material is considered. This material has requires five mechanical constants to characterize its corresponding mechanical behavior. Assuming that each individual ply of the composite is thin enough to be considered as a thin plate/shell element (plane stress state), its matrix of mechanical properties can be described by Eq. (17).

𝑸=

[

𝐸1

𝐸2𝜈12

1 ― 𝜈12𝜈21 𝐸1𝜈12

1 ― 𝜈12𝜈21 𝐸2

1 ― 𝜈12𝜈21

1 ― 𝜈12𝜈21

0

0

0

][

0 𝐺12

]

𝑄11 𝑄12 0 = 𝑄21 𝑄22 0 0 0 𝑄66

(17)

Figure 4: Schematic representation of the fibers orientation relative to the adopted coordinate systems. The (1, 2, 3) system is aligned with the composite fibers for each layer. being 𝐸, 𝐺, and 𝜈 (Poisson ratio) associated to the directions 1 and 2 of Fig. 4. After this procedure, the SHBT approach is applied to determine the equivalent properties to be used in the FE model of the shaft. The advantage of such method is that it permits not only to consider the structural damping and shear stress of the Timoshenko beam elements but also allows for any stacking sequence. This characteristic gives the method a greater field of application, for it does not suffer from the limitation of only representing cross-ply or angleply composites, as in most of the simplified methods. The SHBT approach [24] considers a homogenized flexural stiffness EI with Young’s modulus 𝐸𝑦𝑘 and inertia moment of area 𝐼𝑘 for each ply 𝑘 of the composite material, as shown in Eq. (18).

{

𝜋

𝐼𝑘 = 4(𝑅4𝑘 ― 𝑅4𝑘 ― 1) 𝐸𝑦𝑘 =

[

𝑐𝑜𝑠4(𝜃𝑘) 𝐸1

+

𝑠𝑖𝑛4(𝜃𝑘) 𝐸2

+ 𝑐𝑜𝑠2(𝜃𝑘)𝑠𝑖𝑛2(𝜃𝑘)

(

1 𝐺12

―

)]

2𝜐21 𝐸2

―1

(18)

where Rk is the outer radius of the ply k and Rk-1 is its inner radius. The resulting product 𝐸𝐼 is given by Eq. (19). 𝑁

𝐸𝐼 = ∑𝑘 = 1𝐸𝑦𝑘𝐼𝑘. Equations (20) and (21) present the homogenized shear stress 𝐺𝑆.

(19)

{

𝑆𝑘 = 𝜋(𝑅2𝑘 ― 𝑅2𝑘 ― 1)

[

(

𝐺𝑘12 = 4𝑐𝑜𝑠2(𝜃𝑘)𝑠𝑖𝑛2(𝜃𝑘)

1 𝐸1

1

)

𝜐21

+ 𝐸2 + 2 𝐸1 +

(𝑐𝑜𝑠2(𝜃𝑘) ― 𝑠𝑖𝑛2(𝜃𝑘))2 𝐺𝑘12

]

―1

(20)

𝑁

𝐺𝑆 = ∑𝑘 = 1𝐺𝑘12𝑆𝑘

(21)

Equation (22) is used to obtain the equivalent stiffness that is used in the matrices 𝑫𝑖 and 𝑲𝑖. 𝑁

𝛼𝐸𝐼 = ∑𝑘 = 1𝐸𝑘𝑦𝐼𝑘

(22)

where the damped elastic properties matrix is given by Eq. (23).

{[

[

]

𝜓11𝑄11 𝜓11𝜈21𝑄11 0 𝜓22𝑄22 0 𝑸𝜓 = 𝜓22𝜈21𝑄11 , 𝑖𝑛 (𝑋,𝑌,𝑍), 𝑎𝑛𝑑: 0 0 𝜓12𝑄66

𝑸𝜓 = 𝑻𝑸𝜓𝑻𝑡, 𝑖𝑛 (1, 2, 3), 𝑤ℎ𝑒𝑟𝑒: 𝑐𝑜𝑠 𝜃𝑘) 𝑠𝑖𝑛2(𝜃𝑘) 2𝑠𝑖𝑛(𝜃𝑘)𝑐𝑜𝑠(𝜃𝑘) 2( ) 2( ) (𝜃𝑘)𝑐𝑜𝑠(𝜃𝑘) . 𝑠𝑖𝑛 𝜃 𝑐𝑜𝑠 𝜃 ―2𝑠𝑖𝑛 𝑻= 𝑘 𝑘 2( ) ( ) ( ) ( ) ( ) ―𝑠𝑖𝑛 𝜃𝑘 𝑐𝑜𝑠 𝜃𝑘 𝑠𝑖𝑛 𝜃𝑘 𝑐𝑜𝑠 𝜃𝑘 𝑐𝑜𝑠 𝜃𝑘 ― 𝑠𝑖𝑛2(𝜃𝑘)

]

2(

(23)

𝑘

where ψ is the damping factor. The resulting matrix 𝑸𝜓 is used to determine 𝐸𝑦, as given by Eq. (24).

{

―1

𝑪 = [𝑸𝜓] 𝑘 𝐸𝑦 = [𝑪(2,2)] ―1

(24)

Finally, it can be concluded that all the terms of the equations of motion of a rotor system have been taken into account (see Eq. (2)). Note that the SHBT approach is capable of considering the internal damping effects of the composite material with a good approximation of the influence of the orthotropic layers on the shaft dynamics. The obtained FE model is complemented by the Timoshenko shear stress representation. It is expected that this formulation leads to an improved representation of the dynamic behavior of composite shafts. 3 Rotor test rig Figure 5 presents the rotating machine used in the present contribution. It is composed of a horizontal composite hollow shaft and two aluminum discs (D1 and D2), and supported by two self-aligning ball bearings (B1 and B2). Figure 6a shows a schematic of the test rig together with the main coordinate systems and sensor positions (Prox. H1, Prox. V1, Prox. H2, and Prox. V2). Figure 6b depicts the geometry of the system.

The composite hollow shaft used in the test rig presented in Fig. 5 is made of a special high-modulus and preimpregnated carbon fiber plies, exhibiting twenty layers with the following stacking sequence: [0, 0, 0, 0, 90, 90, 45, -45, 0, 0, 0, 45, -45, 90, 90, 0, 0, 0, 0, 0/90] (degrees of inclination, relative to the Y direction of Fig. 4), from the most internal to the most external ply. The last layer is made of a crossed mesh with interspersed filaments of carbon fibers, where half the fibers are oriented at 0 and the other half at 90 degrees. The geometric properties of the shaft and its plies are shown in Tab. 1. Table 2 presents the physical and geometric properties of the discs. Self-aligning ball bearings were used to support the shaft. They were positioned on the shaft by using sleeves and locking nuts, similar to those used for the discs. The electric motor is driven by a digital controller operated via computer software and was

Figure 5: Composite rotor test rig.

a)

b) Figure 6: Schematic of the experimental test rig for the composite shaft, where: a) configuration of the test rig; b) geometrical configuration and positioning of the discs and bearings.

Table 1: Geometric properties of the shaft. Shaft Properties

Length [m]

Outer diameter [m]

Inner diameter [m]

Density [kg/m³]

Value

0.907

0.018

0.0128

1667

set according to the requirements of the tests performed. This system features a laser encoder effectively measuring the system rotating speed. The rotation of the shaft was set so that no variations greater than ±5 RPM around the operating rotation speeds were observed. Table 2: Physical and geometric properties of the discs. Disks Properties

Value

Thickness [m]

0.016

Outer diameter [m]

0.150

Inner diameter [m]

0.018

Density [kg/m³]

2700

Young’s Modulus [Pa]

69x109

4 Model updating By using the previously described formulation, the SHBT based FE model was assembled to represent the test rig that was discretized by using 39 beam elements. A simple description is given by Fig. 7.

Figure 7: FE model of the experimental test rig. The numbers below the sketch represent the element (with respect to the “left” node) for which the associated component is described. In order to determine the unknown physical properties of the rotor system, two model updating procedures were performed separately. The Differential Evolution (DE) optimization approach was used to solve the associated typical inverse problems (Storn and Price, 1995). The minimization procedure was performed separately due to the number of unknown parameters of the rotating machine. In this case, 19 variables are to be determined for the FE model. The first minimization problem was dedicated to obtaining the shaft parameters, namely Young’s modulus, shear modulus, Poisson’s ratios, and the damping factors ψ (see Eq. (23)). In this case, the chosen objective function was based on the norm of the difference between the experimental and simulated FRFs obtained with the rotor (shaft + discs + ball bearings) supported by wires (i.e., mimicking a free-free condition). The experimental FRFs were

determined by using an electromechanical shaker and two accelerometers (A1 and A2). The shaker and accelerometers were positioned on the shaft along the same direction. Figure 8 presents the experimental setup. The numerical FRFs were obtained by considering the FE model presented in Fig. 7 disregarding the bearings B1 and B2. Table 3 shows the results obtained at the end of the minimization process. Figure 9 depicts the simulated FRFs obtained by using the parameters shown in Tab. 3, and the respective mode shapes for the free-free system configuration. The experimental FRFs are also presented for comparison purposes. Note that the FRFs are similar, demonstrating the representativeness of the obtained shaft parameters. The optimization problem was performed by considering an initial population of 100 individuals and a maximum of 200 iterations.

Figure 8: Composite hollow shaft (shaft + discs + ball bearings) supported by wires. Table 3: Shaft parameters. Variable

Minimum Limit

Optimum Value

Maximum Limit

Young’s modulus 00 [Pa] Young’s modulus 900 [Pa] Young’s modulus 00/900 [Pa] In-plane shear modulus [Pa] In-plane shear modulus 00/900 Major Poisson’s ratio [Pa] Major Poisson’s ratio 00/900 ψ1 ψ2 ψ3

70x109 70x109 20x109 1x108 1x108 0.05 0.05 1x10-8 1x10-8 1x10-8

103.68x109 127.06x109 47.50x109 8.98x108 3.05x109 0.3051 0.2803 2.67x10-7 7.12x10-6 6.58x10-6

150x109 150x109 80x109 10x109 10x109 0.5 0.5 1x10-5 1x10-5 1x10-5

The second minimization problem was dedicated to obtaining the stiffness and damping coefficients of the bearings and the torsional stiffness of the coupling. In this case, the objective function was based on the norm of the difference between the experimental and simulated FRFs obtained for the rotor at rest (shaft + discs + ball bearings) supported by the bearings of the machine, as presented in Fig. 5. The experimental FRFs were determined by using an impact hammer and two proximity probes, namely Prox. H1 and Prox. V1. Impacts were applied separately along the X and Z directions of the discs D1 and D2, respectively. Prox. H1 was used for the impacts along the X direction and Prox. V1 was used for the impacts along the Z direction, resulting 4 FRFs. The numerical FRFs were obtained by considering

the FE model presented in Fig. 7 and the corresponding parameters are shown in Tab. 3. The torsional stiffness of the coupling was applied at node #1 of the FE model (DoFs 𝜃 and 𝜑; see Figs. 1 and 7). In this case, the minimization problem was solved by considering an initial population of 90 individuals and a maximum of 200 iterations for the optimizer. Table 4 shows the results obtained at the end of the minimization process. Figure 10 presents the simulated FRFs obtained by using the parameters shown in Tab. 3 and Tab. 4, and the corresponding mode shapes of the rotating machine. The experimental FRFs are also presented for comparison purposes.

10-2

10-2 Numerical Experimental

Amplitude [m/N]

Amplitude [m/N]

Numerical Experimental

10-4

10-6

10-8

50

100

150

200

10-4

10-6

10-8

250

50

100

Frequency [Hz]

200

a)

b)

ZY plane 1 st mode shape (63.4708 Hz)

ZY plane 2 nd mode shape (180.206 Hz)

250

1

Normalized amplitude ( )

1

Normalized amplitude ( )

150

Frequency [Hz]

0.5

0

-0.5

-1

0.5

0

-0.5

-1 0

0.2

0.4

0.6

Length of the shaft (m)

c)

0.8

1

0

0.2

0.4

0.6

0.8

1

Length of the shaft (m)

d)

Figure 9: Numerical and experimental FRFs obtained for the shaft for a free-free condition, and mode shapes for the first 2 natural frequencies: a) accelerometer A1; b) accelerometer A2; c) 1st mode shape; d) 2nd mode shape. Note that the numerical FRFs obtained along the Z direction are similar to the experimental ones (see Figs. 10b and 10d), particularly in the vicinity of the first natural frequency (26.5 Hz – Z direction). A larger difference can be observed between the numerical and experimental FRFs obtained along the X direction (see Figs. 10a and 10c). In this

case, only the first natural frequency (26.5 Hz – X direction) demonstrated to be satisfactorily identified. Similar results were found by [13]. However, the obtained FE model is considered as being representative enough for the applications found in the context of this work. Figure 11 shows the Campbell diagram of the rotor, in which the first backward whirl speed is, approximately, 1570 RPM, the first forward whirl speed is 1670 RPM, and the second backward whirl speed is 7740 RPM. In this case, only a stable condition was achieved in the considered operating range of the rotating machine (0-9000 RPM; numerical analysis).

10-2

Numerical Experimental

Numerical Experimental

10-4

Amplitude [m/N]

Amplitude [m/N]

10-2

10-6

10-8

10-10

50

100

150

200

10-4

10-6

10-8

10-10

250

50

Frequency [Hz] a)

100

150

Numerical Experimental

Numerical Experimental

10-4

Amplitude [m/N]

Amplitude [m/N]

250

10-2

10-2

10-6

10-8

10-10

200

Frequency [Hz] b)

50

100

150

Frequency [Hz] c)

200

250

10-4

10-6

10-8

10-10

50

100

150

Frequency [Hz] d)

200

250

ZY plane 1 st mode shape (21.0752 Hz)

ZY plane 2 nd mode shape (25.5452 Hz) 1

Normalized amplitude ( )

Normalized amplitude ( )

1

0.5

0

-0.5

-1

0.5

0

-0.5

-1 0

0.2

0.4

0.6

0.8

1

0

0.2

Length of the shaft (m)

0.4

0.6

0.8

1

Length of the shaft (m)

e) f) Figure 10: Numerical and experimental FRFs obtained for the assembled rotor system: a) impacts along the X direction at D2 and Prox. H2; b) impacts along the Z direction at D2 and Prox. V2; c) impacts along the X direction at D1 and Prox. H1; b) impacts along the Z direction at D1 and Prox. V1; e) 1st mode shape; f) 2nd mode shape.

Table 4: Optimized stiffness and damping coefficients. Bearing B1

Bearing B2

Torsional stiffness

Variable

Lower Limit

Optimized Value

Upper Limit

kxx [N/m] kzz [N/m] dxx [Ns/m] dzz [Ns/m] kxx [N/m] kzz [N/m] dxx [Ns/m] dzz [Ns/m] krot [N/rad]

1x102 1x102 1x102 1x102 1x102 1x102 1x102 1x102 1x102

9.9072x104 9.9311x104 2.3510x104 2.2380x104 9.4077x104 9.9883x104 2.3990x104 2.3600x104 2.6713x102

1x107 1x107 1x106 1x106 1x107 1x107 1x106 1x106 1x105

* k: stiffness coefficient; d: damping coefficient.

200

Frequency [Hz]

150

100

50

0

1000 2000 3000 4000 5000 6000 7000 8000 9000

Rotation speed [RPM]

Figure 11: Campbell diagram.

5 Unbalance responses The numerical and experimental unbalance responses of the composite hollow shaft were compared for the rotor system operating at 992 and 1387 RPM. Unbalance masses and their corresponding angular positions at the discs D1 and D2 were identified for both rotation speeds. Then, two optimization processes were applied separately. The chosen objective function was based on the norm of the difference between the experimental and simulated vibration responses of the shaft. The experimental vibration responses were measured by using the proximity probes Prox. H1, Prox. V1, Prox. H2, and Prox. V2. In this case, the minimization problems were solved by considering an initial population of 150 individuals and a maximum of 200 iterations for the optimizer. Table 5 shows the results obtained at the end of the various minimization processes. Figure 12 presents the numerical vibration responses obtained by using the parameters shown in Tab. 5 for the system at 992 RPM. Figure 13 presents the corresponding numerical responses determined for the rotor at 1387 RPM. The experimental vibration responses are also presented for comparison purposes. It is worth mentioning that the vibration responses of the rotating machine were measured by applying order analysis after two hours of operation for each rotation speed. During the experiments, the room temperature was maintained sufficiently stable (set at 20 degrees Celsius). The barely noticeable difference between the numerical and experimental responses for both rotation speeds indicates that Table 5: Unbalance masses and their corresponding angular positions. Variable Unbalance - D1 [Kgm] Angular position - D1 [rad] Unbalance – D2 [Kgm] Angular position - D1 [rad]

Lower Limit Optimized (992 RPM) 1x10-6 0 1x10-6 0

2.10x10-3 4.30 6.72x10-5 3.49

Optimized (1387 RPM)

Upper Limit

1.82x10-3 5.88 4.58x10-5 5.19

1x10-2 2π 1x10-2 2π

the method can determine sufficiently accurate values for the shaft physical properties, independently of the rotating speed considered during the optimization. Thus, the SHBT model is able to properly represent the dynamic behavior of the rotor. Figures 12 and 13 demonstrate that only the 1Ω vibration component of the order signal show relevant amplitude in the spectral responses of the rotor. Thus, the nonlinearity of both composite shaft and damping participate in the responses, however, its influence on the dynamic behavior of the rotor is small at the measured velocities. Additionally, the time domain vibration responses match accordingly, presenting only negligible discrepancies in amplitude and phase for the Z direction. Table 6 presents a comparison between the model and experimental results for the performed analysis.

Figure 12: Comparison of numerical and experimental results both for the order analyses and time signal responses, for the rotation speed of 992 RPM (the sub-indexes indicate the disc and direction, respectively - H stands for Horizontal and V for Vertical directions).

Figure 13: Comparison of numerical and experimental results for the order analysis and time signal response, for the rotation speed of 1387 RPM (the sub-indexes indicate the disc and direction, respectively). Table 6: Comparison of simulated and experimental results: a) phase lag between experimental and numerical time signal; b) simulated amplitude values for the 1Ω peak; c) experimental amplitude values for the 1Ω peak. Disc and direction

a) Phase lag [rad]

b) Numeric amp. [m]

c) Test rig amp. [m]

D1 - Horizontal

0.07

4.002x10-5

4.103x10-5

D1 - Vertical

0.61

6.563x10-5

5.274x10-5

D2 - Horizontal

0.00

1.186x10-4

1.154x10-4

D2 - Vertical

0.00

1.612x10-4

1.450x10-4

D1 - Horizontal

0.05

4.821x10-5

4.837x10-5

D1 - Vertical

0.64

6.392x10-5

6.749x10-5

D2 - Horizontal

0.04

2.537x10-4

2.480x10-4

D2 - Vertical

0.58

2.462x10-4

3.386x10-4

Ω = 992 RPM

Ω = 1387 RPM

Based on the time signals, the simulated and experimental rotor orbits were compared, as shown for each disc and rotation, by keeping the same operating conditions as previously described for the time signal acquisition. Figure 14 presents the corresponding results.

Figure 14: Rotor orbits at both disk positions (numerical and experimental results): a) disk 1 at 992 RPM; b) disk 2 at 992 RPM; c) disk 1 at 1387 RPM; d) disk 2 at 1387 RPM. Notably, the simulated response, in Figs. 12 and 13, tends to be always ahead of the experimental phase response, along with the vertical direction. This effect is believed to be associated with high order directional coupling effects of the composite, which were neglected in the formulation. For the orbits in Fig. 14, this phase lag can be observed, where the simulated orbits are slightly different from the experimental ones. Besides, higher order signals can be seen in the experimental orbits, however, they do not appear in the numerical results. This characteristic has to do with small fabrication issues, such as warping deformation of the shaft, misalignment, residual stress inside the composite mesh, imperfections in the epoxy curing process, bearing clearances, etc. 6 Conclusions In the present contribution important characteristics were introduced into the formulation of the composite hollow shaft: a) a better modeling of the structural damping of the composite material [22], by means of the Kelvin- Voigt representation [28]; b) a more realistic bearing model, accounting for its contribution with localized stiffness and damping, which is an improvement with respect to similar models, such as the [11] overhang rotor; c) possibility of modeling general orientations of the ply fiber [24], without any restriction such as cross-ply or angle-ply; d) capacity of representing the coupling between 𝑋 and 𝑍 directions due to asymmetric bearings and coupling induced by orthotropic properties of the composite. Besides, a formulation that was originally intended for thin-walled (rather than thick-walled composite shafts) was able to satisfactorily predict the dynamic behavior of the experimental test rig. An acceptable agreement between numerical and experimental results was obtained for the tests performed. Even though the proposed methodology accounts for additional physical characteristics of the rotor system, the resulting FE model is still sufficiently simple, so that the use of inverse problem techniques through heuristic

optimization is still affordable. Finally, the computational cost of the proposed methodology is similar to the one found for simplified models, such as the Rayleigh-Ritz and the EMBT models. As further research work, the authors are motivated to investigate new and more complex methodologies to describe the dynamic behavior of thick-walled composite shafts. For this aim, the development of FE models based on the Unified Formulation is scheduled. It is believed that much of the unknown dynamic behavior of composite shafts for rotations above the first critical speed can be better understood through the analysis of the displacement, stress, and strain fields along the thickness direction of each composite layer. 7 Acknowledgments The authors are thankful for the financial support provided to the present research effort by CNPq (574001/20085 and 304546/2018-8), FAPEMIG (TEC-APQ-3076-09, TEC-APQ-02284-15, and TEC-APQ-00464-16), CAPES through the INCT-EIE and PDSE programs, and PETROBRAS. References [1] Arab, S. B.; Rodrigues, J. D.; Bouaziz, S.; Haddar, M. (2017): “A finite element based on Equivalent Single Layer Theory for rotating composite shafts dynamic analysis”, Composite Structures, Volume 178, Pages 135144. [2] Arab, S. B.; Rodrigues, J. D.; Bouaziz, S.; Haddar, M. (2017): “Dynamic analysis of laminated rotors using a layerwise theory”, Composite Structures, Volume 182, Pages 335-345. [3] Barbosa, P. C. P. F.; Cavalini Jr., A. A.; Steffen Jr., V. (2018): “Analysis of the dynamic behavior of composite shafts on rotating machines”. 43p. Master’s thesis, Federal University of Uberlândia, Brazil. [4] Brush, M. (1999): “Still Spinning After All These Years: A Profile of the Ultracentrifuge”. The Scientist, v. 13, p. 16-18, Out. [5] Bucciarelli, L. L. (1982): “On the Instability of Rotating Shafts Due to Internal Damping”, Journal of Applied Mechanics, Volume 49, Issue 2, Pages 425-428. [6] Cavalini Jr., A.; Guimarães, T.; Da Silva, B. R. M. G.; Steffen Jr., V. (2017): “Analysis of the Dynamic Behavior of a Rotating Composite Hollow Shaft”, Latin American Journal of Solids and Structures, Volume 14, Pages 1-16. [7] Cavalini Jr, A. A. (2013): “Detection and identification of incipient transversal cracks in flexible and horizontal shafts of rotating machines”. 270p. PhD Thesis, Federal University of Uberlândia, Brazil. [8] Del Claro, V. T. S.; Barbosa, P. C. P. F.; Cavalini Jr., A. A.; Steffen Jr., V. (2017): “A shell based fem model for thick walled composite rotors”, COBEM-2017-2118, 24th ABCM International Congress of Mechanical Engineering, December 3-8, Curitiba, PR, Brazil. [9] Duchemin M. (2003): “Contribution à l’étude du comportement dynamique d’un rotor embarqué”. PhD thesis, INSA Lyon, Lyon. [10] Duchemin, M.; Berlioz, A.; Ferraris, G. (2006): “Dynamic behavior and stability of a rotor under base excitation”. ASME J Vib Acoust 128:576.

[11] Friswell, M. I.; Penny, E. T.; Garvey, D.; Lees, W. (2010): “Dynamics of Rotating Machines”, Cambridge Aerospace Series, Cambridge University Press, 544p. ISBN: 978-0-521-85016-2. [12] Gubran, H. B. H.; Gupta, K. (2005): “The effect of stacking sequence and coupling mechanisms on the natural frequencies of composite shafts”, Journal of Sound and Vibrations, Volume 282(1-2), Pages 231-248. [13] Gupta, K. (2015): “Composite Shaft Rotor Dynamics: An Overview”. In: Proceedings of VETOMAC, X, Manchester. Vibration Engineering and Technology of Machinery. p. 79-94. [14] Lalanne, M.; Ferraris, G. (1998): “Rotordynamics Prediction in Engineering”. New York: J. Wiley and Sons, 266p. [15] Lobato, F. S.; Steffen Jr., V. (2017): “Multi-Objective Optimization Problems: Concepts and Self-Adaptive Parameters with Mathematical and Engineering Applications”, Springer International Publishing, 160p. ISBN: 978-3-319-58565-9. [16] Mendonça, W. R. P.; Medeiros, E. C.; Pereira, A. L. R.; Mathias, M. H. (2017): ” The dynamic analysis of rotors mounted on composite shafts with internal damping”, Composite Structures, Volume 167, Pages 50-62. [17] Milazzo, A. (2014): “Refined equivalent single layer formulations and finite elements for smart laminates free vibrations”, Composites Part B: Engineering, Volume 61, Pages 238-253. [18] Milazzo, A. (2014): “Layer-wise and equivalent single layer models for smart multilayered plates”, Composites Part B: Engineering, Volume 67, Pages 62-75. [19] Reddy, J. N. (2017): “Energy Principles and Variational Methods in Applied Mechanics”, Wiley, 3rd Edition, 760p. ISBN: 978-1-119-08737-3. [20] Silveira, M. E. (2001): “Análise do Comportamento Dinâmico de Rotores em Eixos Bobinados”. 97 p. Master’s thesis, Federal University of Santa Catarina, Florianópolis, Brazil. [21] Singh, S.P.; Gupta, K. (1996): “Composite shaft rotordynamic analysis using a layerwise theory”, Journal of Sound and Vibration, Volume 191(5), Pages 739-756. [22] Singh, S. P.; Gupta, K. (1994): “Free damped flexural vibration analysis of composite cylindrical tubes using beam and shell theories”, Journal of Sound and Vibration, Volume 172(2), Pages 171-190. [23] Sino, R.; Baranger, A. B.; Chatelet, A. G.; Jacquet, A. (2008): “Dynamic analysis of a rotating composite shaft”, Composites Science and Technology, Volume 68, Pages 337–345. [24] Sino, R. (2007): “Comportement Dynamique et Stabilité des Rotors: Application aux Rotors Composites”. 157 p. PhD Thesis, INSA - Lyon, Lyon. [25] Sousa Jr., M. S.; Cavalini Jr., A. A.; Steffen Jr., V. (2017): “Analysis of the Dynamic Behavior of Onboard Rotor”. 50 p. Master’s thesis, Federal University of Uberlândia, Brazil. [26] Sousa Jr., M. S.; Del Claro, V. T. S.; Cavalini Jr., A. A.; Steffen Jr., V. (2017):” Numerical investigation on the dynamic behavior of an onboard rotor system by using the FEM approach”, Journal of the Brazilian Society for Mechanical Sciences and Engineering, Volume 39, Pages 2447-2458. [27] Tsai, S. W. (1988): “Composites design”, 4th Edition, Ohio, USA, Dayton. [28] Wasilkoski, C. M. (2006): “Comportamento Mecânico dos Materiais Poliméricos”. 82p. PhD Thesis, Federal University of Paraná, Brazil.

[29] Wettergreen, H.L.; Olsson, K.O. (1996): “Dynamic Instability of a Rotating Asymmetric Shaft with Internal Viscous Damping Supported in Anisotropic Bearings”. Journal of Sound and Vibration, v. 195, p. 75-84, Ago. [30] Zhao, J.; Wong, P. K.; Ma, X.; Xie, Z.; Xu, J.; Cristino, V.A. (2018):"Simplification of finite element modeling for plates structures with constrained layer damping by using single-layer equivalent material properties", Composites Part B: Engineering.