Dynamic modeling and coupling characteristics of rotating inclined beams with twisted-shape sections

Abstract

In the existing literature, most studies investigated the free vibrations of a rotating pre-twisted cantilever beam; however, few considered the effect of the elastic-support boundary and the quantification of modal coupling degree among different vibration directions. In addition, Coriolis, spin softening, and centrifugal stiffening effects are not fully included in the derived equations of motion of a rotating beam in most literature, especially the centrifugal stiffening effect in torsional direction. Considering these deficiencies, this study established a coupled flapwise–chordwise–axial–torsional dynamic model of a rotating double-tapered, pre-twisted, and inclined Timoshenko beam with elastic supports based on the semi-analytic method. Then, the proposed model was verified with experiments and ANSYS models using Beam188 and Shell181 elements. Finally, the effects of setting and pre-twisted angles on the degree of coupling among flapwise, chordwise, and torsional directions were quantified via modal strain energy ratios. Results showed that 1) the appearance of torsional vibration originates from the combined effect of flapwise–torsional and chordwise–torsional couplings dependent on the Coriolis effect, and that 2) the flapwise–chordwise coupling caused by the pure pre-twisted angle is stronger than that caused by the pure setting angle.

Key words： elastic-support boundary; pre-twisted beam; semi-analytic method; modal strain energy ratio; torsional vibration

Introduction

Rotating pre-twisted beam-like structures are extensively used in machinery; aerospace industry; and power systems [¹], such as fan blades [²], marine propellers [³], and wind turbines [⁴]. The important feature of this kind of structure is that its modal characteristics are closely related to many engineering failure problems such as flutter, stability, and fatigue failure [⁵,⁶]. Therefore, many studies have investigated its dynamic characteristics through experimental, analytical, and numerical methods [⁷–¹²].

Early studies were primarily focused on the rotating pre-twisted cantilever beam. Adopting Hamilton’s principle and Newtonian method, Hodges and Dowell [¹³] theoretically derived the coupled flapwise–chordwise–axial–torsional equations of motion of a rotating pre-twisted Euler–Bernoulli beam and found that Hamilton’s principle was easier and more precise to handle the equations of motion than the Newtonian method due to lesser chance of inadvertently omitting some important terms. Zhu [¹⁴] utilized the Rayleigh–Ritz method to establish the coupled flapwise–chordwise–axial equations of motion of a pre-twisted rotating Timoshenko beam and discussed the effects of some parameters (i.e., slenderness ratio, hub radius ratio, and rotational speed) on the natural frequencies of the structure. Ma et al. [¹⁵] discussed the effect of the tip shroud rubbing on the coupled flapwise–axial vibration responses of the rotating untwisted Euler–Bernoulli beam. In their later work [¹⁶], the effects of the pre-twisted angle, setting angle, and shear deformation were further included. However, the developed model in Ref. [¹⁶] was not suitable for simulating the dynamic behaviors of the structure with a larger pre-twisted angle. Sinha [¹⁷] utilized the Rayleigh–Ritz method to determine the coupled flapwise–axial–torsional vibration characteristics of a rotating pre-twisted beam suffering from contact load and Coulomb friction. However, the limitation of the proposed model in Ref. [¹⁷] was its inability to acquire the classical two-stripe mode of the structure. Yang et al. [¹⁸] adopted the power series method to investigate the dynamic frequencies varying with the rotating speed and setting angles, as well as corresponding complex mode shapes of a rotating tapered and pre-twisted cantilever beam with flapwise–chordwise–axial–torsional coupling. Şakar and Sabuncu [¹⁹] used a finite element method to establish the coupled flapwise–chordwise–axial dynamic model of a rotating pre-twisted aerofoil cross-section cantilever Euler–Bernoulli beam. However, both spin softening and Coriolis effects were ignored in Ref. [¹⁹]. Based on the same mathematical model in Ref. [¹⁹], the effect of shear deformation was further considered by Sabuncu and Evran [²⁰]. They found that the stability of the studied structure increased with decreasing stagger angle and increasing rotational speed and disk radius. Subrahmanyam et al. [²¹] utilized the Galerkin method in combination with a linear perturbation technique to solve the coupled flapwise–chordwise–torsional equations of motion of a rotating pre-twisted Euler–Bernoulli cantilever beam. The authors claimed that the accuracy of the structure with the second-degree geometric nonlinearities included was adequate and the introduction of nonlinear terms into the thin blades easily resulted in torsional divergence. Discarding the effects of Coriolis and axial inertia, Sina and Haddadpour [²²] established the torsional equations of motion of a rotating pre-twisted thin-walled composite beam and indicated that the pre-twisted angle and material anisotropy had remarkable influence on the torsional behaviors of the studied structure, including torsional frequency and hardening/softening effect. Adair and Jaeger [²³] ignored the spin softening and Coriolis effects and then applied the modified Adomian decomposition method to obtain the coupled flapwise–chordwise free vibration characteristics of a uniform pre-twisted rotating Euler–Bernoulli beam. The developed model was verified using natural frequency comparisons with the published results, and the influence of the pre-twisted angle and rotating speed on the frequency characteristics was also investigated. Oh and Yoo [²⁴] combined Kane’s method with the Galerkin method to establish the coupled flapwise–chordwise–axial–torsional dynamic model of a rotating pre-twisted cantilever beam with arbitrary cross-section. The numerical results indicated that the proposed model had good precision via the comparisons of the modal characteristics obtained from a commercial finite element code. Lee and Lee [²⁵] and Banerjee [²⁶] applied the transfer matrix method to determine the exact flapwise and chordwise modal characteristics of a rotating pre-twisted cantilever beam, and the availability of the proposed models was also verified using natural characteristic comparisons with the literatures. Subrahmanyam and Kaza [²⁷] employed the finite-difference approach and the potential energy method to solve the coupled bending–bending equations of motion of a torsionally rigid slender beam. The results showed that the inclusion of the Coriolis effects was significant for beams with moderate to large thickness ratios and insignificant for beams with small thickness ratios. Hashemi and Richard [²⁸] utilized a dynamic finite element to determine the natural frequencies and mode shapes of rotating assemblages made of beams, which was verified via two illustrative examples of vertical and radial beams with Coriolis effect included. Banerjee and Kennedy [²⁹] used an exact dynamic stiffness method to study the in-plane free vibrations of a rotating beam, which covered the effects of Coriolis force, hub radius, and outboard force. The results revealed that the coupling between bending and axial deformations was rather distinct. Considering the effects of anisotropic material, rotary inertia, and shear deformation, Oh et al. [³⁰] established the coupled bending–bending equations of motion of a rotating pre-twisted box beam, which was verified by means of eigenfrequency characteristics against the results obtained from the literatures. Latalski et al. [³¹] investigated the coupled bending–shear–twist vibration characteristics of a thin-walled rotating cantilevered beam composed of composite material. In particular, the hub equation resulting from the rotation angle was also included in the structure. The results showed that the maximum magnitude of the flexural–torsional mode coupling for the circumferentially asymmetric stiffness configuration occurred at a ply angle of 74°. Ondra and Titurus [³²] established the coupled bending–bending–torsional equations of motion of the rotating pre-twisted beam–tendon system, whose natural frequencies and mode shapes are determined via the combination of a boundary value problem solver and differential quadrature method. In addition, the effects of rotation, pre-twist, and cross-sectional coupling on the modal characteristics of the system were also discussed.

However, in many engineering applications, it is very difficult to enforce the infinite rigidity at the cantilever end of the pre-twisted beam-like structure. In Refs. [³³–³⁵], the numerical and experimental results indicated some differences in the low-order modes between a pre-twisted cantilever blade with a tenon and that without a tenon, which also indicated that the root flexibility for the practical cantilevered beam-like structures should be considered. Ignoring the Coriolis effect, the coupled flapwise–chordwise dynamic characteristics of a rotating, in-extensional, and pre-twisted Euler–Bernoulli beam with elastic free boundaries were solved by Lin [³⁶] using the transition matrix method. Based on the developed model, the instability mechanism of the rotating beam was investigated under different parameter influences. In their later works [³⁷,³⁸], a tip mass was added into the preceding model in Ref. [³⁶]. In Ref. [³⁷], the self-adjointness of the structure and the symmetric properties of the matrix Green functions were revealed, and some parameter influences on the natural frequencies of pre-twisted and untwisted beams are discussed in Ref. [³⁸]. Ignoring the Coriolis effect, Choi and Chou [³⁹] applied the modified differential quadrature method to solve the coupled flapwise–chordwise free vibration of a rotating, in-extensional, and pre-twisted Timoshenko beam with elastic–elastic boundaries. The advantages of the proposed method were the permission of any combinations of end restraining conditions. Bambill et al. [⁴⁰] utilized the differential quadrature method in combination with the domain decomposition to analyze the flapwise vibration behaviors of a rotating tapered multi-span beam under the effect of spring stiffness. The results showed that the translational spring stiffness had more significant effects on the dynamic characteristics of the structure than the rotational spring stiffness. Digilov and Abramovich [⁴¹] conducted theoretical and experimental studies of the impact of root flexibility on the flapwise vibration of a uniform beam where the shear deformation and rotary inertia were left out.

According to the literatures listed above, most were conducted on the free vibrations of the rotating pre-twisted beam with clamped boundary conditions, whose torsional vibration are ignored in most cases. In addition, the centrifugal stiffening effect in the torsional direction and Coriolis effect were overlooked in many literatures. In particular, the introduction of non-zero setting angle, pre-twisted angle, and rotating speed can result in vibration coupling among the flapwise, chordwise, and torsional directions, but quantitative studies on the corresponding degree of coupling are limited. Considering these deficiencies, the present study established a dynamic model of the rotating inclined beams with twisted-shape sections and elastically restrained root using Hamilton’s principle in combination with the assumed mode method. Numerical simulation in ANSYS and experimental test were conducted to verify the proposed model. Next, the modal strain energy ratio (MSER) affected by the setting angle, pre-twisted angle, and rotating speed was utilized to quantify the degree of coupling among flapwise, chordwise, and torsional directions. Finally, some conclusions were made.

Mathematical modeling

Establishment of kinetic and potential energies

In engineering applications, it is difficult to apply the infinite rigidity to the beam root, thereby leading to the existence of certain flexibilities. On this basis, the schematic of a pre-twisted and double-tapered beam with elastic supports is shown in Fig. 1. In Fig. 1, OXYZ, O^rX^rY^rZ^r, o^dx^dy^dz^d, and oxyz are the global coordinate system, the rotating coordinate system, the local coordinate system attached to the joint surface on the disk, and the local coordinate system attached to the arbitrary beam section, respectively; β₀ and γ(L) are the initial setting and pre-twisted angles, respectively (see Fig. 1(b)); three linear (u, v, w) and three angular (θ, φ, ψ) displacement components of o in oxyz are assumed (see Fig. 1(c)). Accordingly, three linear (k_x, k_y, k_z) and three angular (k_rx, k_ry, k_rz) supporting stiffness are established to simulate the flexibility of the root section of the beam, as shown in Fig. 1(a).

Fig.1 Pre-twisted and double-tapered beam: (a) Reference frames, (b) pre-twisted and setting angles, and (c) deformation modes of an arbitrary section.

Full size|PPT slide

In Fig. 1(b), the setting angle β of an arbitrary beam section can be calculated based on the linear interpolation hypothesis:

(1)\beta ={{\displaystyle \beta }}_{\mathrm{0}}+\mathit{\gamma },

where γ=γ(L)x/L is the pre-twisted angle at the initial coordinate x in o^dx^dy^dz^d and L is the beam length.

Combining Fig. 1(a) with Fig. 1(c), the coordinates of an arbitrary point “p” on the beam section in OXYZ can be expressed as follows:

(2){\mathit{r}}_p=\left[\begin{array}{c}\hfill X_p\hfill \\ \hfill Y_p\hfill \\ \hfill Z_p\hfill \end{array}\right]={\mathit{A}}_{\mathrm{5}}{\mathit{A}}_{\mathrm{4}}\left(\left[\begin{array}{c}\hfill R_{\mathrm{d}}+x+u\hfill \\ \hfill v\hfill \\ \hfill w\hfill \end{array}\right]+{\mathit{A}}_{\mathrm{3}}{\mathit{A}}_{\mathrm{2}}{\mathit{A}}_{\mathrm{1}}\left(\begin{array}{c}\hfill \mathrm{0}\hfill \\ \hfill y\hfill \\ \hfill z\hfill \end{array}\right)\right),

where R_d is the disk radius, y and z are the coordinates of p in oyz, and A₁, A₂, A₃, A₄, and A₅ are a series of transformation matrices from oxyz to OXYZ, whose expressions are shown as follows:

{\mathit{A}}_{\mathrm{1}}=\left[\begin{array}{ccc}\hfill \mathrm{1}\hfill & \hfill -\psi \hfill & \hfill \mathrm{0}\hfill \\ \hfill \psi \hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill \end{array}\right],{\mathit{A}}_{\mathrm{2}}=\left[\begin{array}{ccc}\hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \varphi \hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill \\ \hfill -\varphi \hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill \end{array}\right],

{\mathit{A}}_{\mathrm{3}}=\left[\begin{array}{ccc}\hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill -\theta \hfill \\ \hfill \mathrm{0}\hfill & \hfill \theta \hfill & \hfill \mathrm{1}\hfill \end{array}\right],{\mathit{A}}_{\mathrm{4}}=\left[\begin{array}{ccc}\hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \cos \beta \hfill & \hfill -\sin \beta \hfill \\ \hfill \mathrm{0}\hfill & \hfill \sin \beta \hfill & \hfill \cos \beta \hfill \end{array}\right],

(3){\mathit{A}}_{\mathrm{5}}=\left[\begin{array}{ccc}\hfill \cos \alpha \hfill & \hfill -\sin \alpha \hfill & \hfill \mathrm{0}\hfill \\ \hfill \sin \alpha \hfill & \hfill \cos \alpha \hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill \end{array}\right].

Here, the small deformation in this work is assumed, and only first-order Taylor expansion is adopted in Eq. (3).

To make the mode shapes of the uniform straight link/beam applicable to discretize the equations of motion of the pre-twisted beam, some necessary displacement transformations are needed first. In this work, the displacements q = [u, v, w, θ, φ, ψ]^T in oxyz are represented by the displacement Q = [U, V, W, Q, F, Y]^T in o^dx^dy^dz^d.

(4)\mathit{q}=\left[\begin{array}{cccccc}\hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \cos \mathit{\gamma }\hfill & \hfill \sin \mathit{\gamma }\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill -\sin \mathit{\gamma }\hfill & \hfill \cos \mathit{\gamma }\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \cos \mathit{\gamma }\hfill & \hfill \sin \mathit{\gamma }\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill -\sin \mathit{\gamma }\hfill & \hfill \cos \mathit{\gamma }\hfill \end{array}\right]\mathit{Q}.

Substituting Eqs. (3) and (4) into Eq. (2) the kinetic energy T_b of the rotating, pre-twisted, double-tapered, and elastic-supported beam can be calculated as follows:

(5)T_{\mathrm{b}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\rho {\displaystyle {\int }_{\mathrm{0}}^L\mathrm{d}x{\displaystyle \underset{A}{\int \int }{\dot{\mathit{r}}}_p^{\mathrm{T}}}}{\dot{\mathit{r}}}_p\mathrm{d}y\mathrm{d}z,

where an over-dot represents the differentiation with respect to time t, r is the density, and A is the sectional area at x in o^dx^dy^dz^d (see Fig. 1(a)).

The corresponding potential energy U_b is given as follows:

(6)U_{\mathrm{b}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}E{\int }_{\mathrm{0}}^L[{AU}^{\prime \mathrm{2}}+(I_z{\cos }^{\mathrm{2}}\mathit{\gamma }+I_y{\sin }^{\mathrm{2}}\mathit{\gamma }){\Psi }^{\prime \mathrm{2}}+(I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right){\Psi }^{\prime }{\mathrm{\Phi }}^{\prime }+(I_z\sin {}^{\mathrm{2}}\mathit{\gamma }+I_y{\cos }^{\mathrm{2}}\mathit{\gamma }){\mathrm{\Phi }}^{\prime \mathrm{2}}+{\displaystyle \frac{GJ}{E}}{\mathrm{\Theta }}^{\prime \mathrm{2}}]\mathrm{d}x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\rho {\dot{\alpha }}^{\mathrm{2}}{\int }_{\mathrm{0}}^{L}A(R_{\mathrm{d}}+x)\left[{\int }_{\mathrm{0}}^x(W^{\prime \mathrm{2}}+V^{\prime \mathrm{2}}+\hat{{\displaystyle \frac{I_x}{A}}{\mathrm{\Theta }}^{\prime \mathrm{2}}})\mathrm{d}x\right]\mathrm{d}x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}G{\int }_{\mathrm{0}}^L{\kappa }_{y}A{(V^{\prime }-\Psi )}^{\mathrm{2}}\mathrm{d}x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}G{\int }_{\mathrm{0}}^L{\kappa }_{z}A{(W^{\prime }+\mathrm{\Phi })}^{\mathrm{2}}\mathrm{d}x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{k}_x{U}^{\mathrm{2}}{|}_{x=\mathrm{0}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{k}_y{V}^{\mathrm{2}}{|}_{x=\mathrm{0}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{k}_z{W}^{\mathrm{2}}{|}_{x=\mathrm{0}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{k}_{\mathrm{r}x}{\mathrm{\Theta }}^{\mathrm{2}}{|}_{x=\mathrm{0}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{k}_{\mathrm{r}y}{\mathrm{\Phi }}^{\mathrm{2}}{|}_{x=\mathrm{0}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{k}_{\mathrm{r}z}{\Psi }^{\mathrm{2}}{|}_{x=\mathrm{0}},

where E and G are Young’s modulus and shear modulus, respectively, I_y=h{b}^{\mathrm{3}}/\mathrm{12}and I_z=b{h}^{\mathrm{3}}/\mathrm{12} are the second moments of area in the y- and z-axes at x in o^dx^dy^dz^d, respectively, b=\left(\mathrm{1}-{\tau }_{\mathrm{b}}x/L\right)b_{\mathrm{0}} and h=\left(\mathrm{1}-{\tau }_{\mathrm{h}}x/L\right)h_{\mathrm{0}} are the width and thickness of the beam section located at x in o^dx^dy^dz^d, b₀ and h₀ are the width and thickness of the root section (see Fig. 1(a)), respectively, t_b and t_h are the tapered ratios of the width and thickness, respectively; A = bh is the area of the beam section at x in o^dx^dy^dz^d, and k_y and k_z are the shear coefficients in the y- and z-directions, respectively. In Eq. (6), \widehat{(·)} are the terms related to centrifugal stiffening effect in torsional vibration [⁴²].

In Eq. (6), J is the torsional moment of inertia, and the corresponding computational formula is written as follows [⁴³]:

J=J_{\mathrm{s}}+J_{\mathrm{a}},J_{\mathrm{s}}=\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}-{\displaystyle \frac{\mathrm{64}h}{{\mathrm{\pi }}^{\mathrm{5}}b}}\tanh \left({\displaystyle \frac{\mathrm{\pi }b}{\mathrm{2}h}}\right)\right]b{h}^{\mathrm{3}},

(7)J_a={\displaystyle \frac{E\mathit{\gamma }{\left(L\right)}^{\mathrm{2}}}{\mathrm{720}G{L}^{\mathrm{2}}}}\left[bh\left(\mathrm{9}{b}^{\mathrm{4}}+\mathrm{10}{b}^{\mathrm{2}}{h}^{\mathrm{2}}+\mathrm{9}{h}^{\mathrm{4}}\right)-\mathrm{60}{I}_x\left(b^{\mathrm{2}}+h^{\mathrm{2}}\right)\right],

where b and h are the width and thickness of an arbitrary rectangular beam section, respectively, and I_x = I_y + I_z is the polar moment of inertia.

Galerkin discretization

Applying Hamilton’s principle to Eqs. (5) and (6), six partial differential equations of motion can be then determined as follows:

Partial differential equations of motion in the x^d-direction:

(8)\left\{\begin{array}{l}{\displaystyle {\int }_{\mathrm{0}}^L\rho A}\ddot{U}\mathrm{\delta }U\mathrm{d}x-{\displaystyle {\int }_{\mathrm{0}}^L\rho A{\dot{\alpha }}^{\mathrm{2}}U\mathrm{\delta }U\mathrm{d}x}\hfill \\ -\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho A\dot{\alpha }\cos {\beta }_{\mathrm{0}}\dot{V}\mathrm{\delta }U\mathrm{d}x}+\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho A\sin {\beta }_{\mathrm{0}}\dot{\alpha }\dot{W}}\mathrm{\delta }U\mathrm{d}x\hfill \\ -{\displaystyle {\int }_{\mathrm{0}}^L\rho A\ddot{\alpha }\cos {\beta }_{\mathrm{0}}V\mathrm{\delta }U\mathrm{d}x}+{\displaystyle {\int }_{\mathrm{0}}^L\rho A\ddot{\alpha }\sin {\beta }_{\mathrm{0}}}W\mathrm{\delta }U\mathrm{d}x\hfill \\ +E{\displaystyle {\int }_{\mathrm{0}}^{L}A{U}^{\prime }\mathrm{\delta }}{U}^{\prime }\mathrm{d}x-{\displaystyle {\int }_{\mathrm{0}}^L\rho A\left(R_{\mathrm{d}}+x\right){\dot{\alpha }}^{\mathrm{2}}\mathrm{\delta }U\mathrm{d}x+k_{x}U\mathrm{\delta }U|{}_{x=\mathrm{0}}}\hfill \end{array}\right\}=\mathrm{0.}

Partial differential equations of motion in the y^d-direction:

(9)\left\{\begin{array}{l}{\displaystyle {\int }_{\mathrm{0}}^L\rho A\ddot{V}}\mathrm{\delta }V\mathrm{d}x-{\displaystyle {\int }_{\mathrm{0}}^L\rho A{\dot{\alpha }}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\beta }_{\mathrm{0}}}V\mathrm{\delta }V\mathrm{d}x\hfill \\ +{\displaystyle {\int }_{\mathrm{0}}^L\rho A\ddot{\alpha }\cos {\beta }_{\mathrm{0}}U\mathrm{\delta }V}\mathrm{d}x\hfill \\ +\mathrm{0.5}{\displaystyle {\int }_{\mathrm{0}}^L\rho A{\dot{\alpha }}^{\mathrm{2}}\sin \left(\mathrm{2}{\beta }_{\mathrm{0}}\right)W\mathrm{\delta }V\mathrm{d}x}\hfill \\ +\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho A\dot{\alpha }\cos {\beta }_{\mathrm{0}}\dot{U}\mathrm{\delta }V\mathrm{d}x}\hfill \\ +{\displaystyle {\int }_{\mathrm{0}}^L\rho A\ddot{\alpha }\cos {\beta }_{\mathrm{0}}\left(R_{\mathrm{d}}+x\right)\mathrm{\delta }V\mathrm{d}x}\hfill \\ +\rho {\dot{\alpha }}^{\mathrm{2}}{\displaystyle {\int }_{\mathrm{0}}^L\left(R_{\mathrm{d}}+x\right)A\left({\displaystyle {\int }_{\mathrm{0}}^x{V}^{\prime }\mathrm{\delta }{V}^{\prime }\mathrm{d}x}\right)\mathrm{d}x}\hfill \\ +{\kappa }_{y}G{\displaystyle {\int }_{\mathrm{0}}^{L}A\left(V^{\prime }-\Psi \right)\mathrm{\delta }{V}^{\prime }\mathrm{d}x}+k_{y}V\mathrm{\delta }V|{}_{x=\mathrm{0}}\hfill \end{array}\right\}=\mathrm{0.}

Partial differential equations of motion in the z^d-direction:

(10)\left\{\begin{array}{l}{\displaystyle {\int }_{\mathrm{0}}^L\rho A}\ddot{W}\mathrm{\delta }W\mathrm{d}x-{\displaystyle {\int }_{\mathrm{0}}^L\rho A{\dot{\alpha }}^{\mathrm{2}}{\sin }^{\mathrm{2}}{\beta }_{\mathrm{0}}W\mathrm{\delta }W\mathrm{d}x}\hfill \\ -{\displaystyle {\int }_{\mathrm{0}}^L\rho A\ddot{\alpha }\sin {\beta }_{\mathrm{0}}U\mathrm{\delta }W\mathrm{d}x}+\mathrm{0.5}{\displaystyle {\int }_{\mathrm{0}}^L\rho A{\dot{\alpha }}^{\mathrm{2}}\sin \left(\mathrm{2}{\beta }_{\mathrm{0}}\right)V\mathrm{\delta }W\mathrm{d}x}\hfill \\ -\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho A\dot{\alpha }\sin {\beta }_{\mathrm{0}}\dot{U}\mathrm{\delta }W\mathrm{d}x}-{\displaystyle {\int }_{\mathrm{0}}^L\rho A\ddot{\alpha }\left(R_{\mathrm{d}}+x\right)\sin {\beta }_{\mathrm{0}}\mathrm{\delta }W\mathrm{d}x}\hfill \\ +{\displaystyle {\int }_{\mathrm{0}}^L\rho {\dot{\alpha }}^{\mathrm{2}}A\left(R_{\mathrm{d}}+x\right)\left({\displaystyle {\int }_{\mathrm{0}}^x{W}^{\prime }\mathrm{\delta }{W}^{\prime }\mathrm{d}x}\right)\mathrm{d}x}\hfill \\ +k_{z}G{\displaystyle {\int }_{\mathrm{0}}^{L}A\left(W^{\prime }+\mathrm{\Phi }\right)\mathrm{\delta }{W}^{\prime }\mathrm{d}x}+k_{z}W\mathrm{\delta }W{|}_{x=\mathrm{0}}\hfill \end{array}\right\}=\mathrm{0.}

Partial differential equations of motion in the rotx^d-direction:

(11)\left\{\begin{array}{l}{\int }_{\mathrm{0}}^L\rho (I_z+I_y)\ddot{\mathrm{\Theta }}\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x-{\int }_{\mathrm{0}}^L\rho {\dot{\alpha }}^{\mathrm{2}}(I_y{\cos }^{\mathrm{2}}\beta +I_z{\sin }^{\mathrm{2}}\beta )\mathrm{\Theta }\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x\hfill \\ -\mathrm{2}{\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \mathit{\gamma }\cos \beta +I_z\sin \mathit{\gamma }\sin \beta )\dot{\mathrm{\Phi }}\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x\hfill \\ -\mathrm{2}{\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \beta \sin \mathit{\gamma }-I_z\cos \mathit{\gamma }\sin \beta )\dot{\Psi }\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x\hfill \\ -{\int }_{\mathrm{0}}^L\rho \ddot{\alpha }(I_z\sin \mathit{\gamma }\sin \beta +I_y\cos \mathit{\gamma }\sin \beta )\mathrm{\Theta }\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x\hfill \\ +{\int }_{\mathrm{0}}^L\rho \ddot{\alpha }(I_z\cos \mathit{\gamma }\sin \beta -I_y\sin \mathit{\gamma }\cos \beta )\Psi \mathrm{\delta }\mathrm{\Theta }\mathrm{d}x\hfill \\ +\mathrm{0.5}{\int }_{\mathrm{0}}^L\rho {\dot{\alpha }}^{\mathrm{2}}(I_z-I_y)\sin \left(\mathrm{2}\beta \right)\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x+G{\int }_{\mathrm{0}}^{L}J{\mathrm{\Theta }}^{\prime }\mathrm{\delta }{\mathrm{\Theta }}^{\prime }\mathrm{d}x\hfill \\ +\rho {\dot{\alpha }}^{\mathrm{2}}{\int }_{\mathrm{0}}^L(R_{\mathrm{d}}+x)A\left({\int }_{\mathrm{0}}^x{\displaystyle \frac{I_x}{A}}{\mathrm{\Theta }}^{\prime }\mathrm{\delta }{\mathrm{\Theta }}^{\prime }\mathrm{d}x\right)\mathrm{d}x+k_{\mathrm{r}\mathit{x}}\mathrm{\Theta }\mathrm{\delta }\mathrm{\Theta }{|}_{x=\mathrm{0}}\hfill \end{array}\right\}=\mathrm{0}.

Partial differential equations of motion in the roty^d-direction:

(12)\left\{\begin{array}{l}{\int }_{\mathrm{0}}^L\rho (I_y{\cos }^{\mathrm{2}}\mathit{\gamma }+I_z{\sin }^{\mathrm{2}}\mathit{\gamma })\ddot{\mathrm{\Theta }}\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x\hfill \\ -{\int }_{\mathrm{0}}^L\rho {\dot{\alpha }}^{\mathrm{2}}(I_y{\cos }^{\mathrm{2}}\mathit{\gamma }+I_z{\sin }^{\mathrm{2}}\mathit{\gamma })\mathrm{\Phi }\mathrm{\delta }\mathrm{\Phi }\mathrm{d}x\hfill \\ +\mathrm{2}{\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \mathit{\gamma }\cos \beta +I_z\sin \mathit{\gamma }\sin \beta )\dot{\mathrm{\Theta }}\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x\hfill \\ +{\int }_{\mathrm{0}}^L\rho \ddot{\alpha }(I_y\cos \mathit{\gamma }\cos \beta +I_z\sin \mathit{\gamma }\sin \beta )\mathrm{\Theta }\mathrm{\delta }\mathrm{\Phi }\mathrm{d}x\hfill \\ +{\mathrm{0.5}\int }_{\mathrm{0}}^L\rho (I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right)\ddot{\Psi }\mathrm{\delta }\mathrm{\Phi }\mathrm{d}x\hfill \\ +{\int }_{\mathrm{0}}^L\rho [I_y\cos \mathit{\gamma }\sin \beta -I_z\sin \mathit{\gamma }\cos \beta ]\ddot{\alpha }\mathrm{\delta }\mathrm{\Phi }\mathrm{d}x\hfill \\ +\mathrm{0.5}E{\int }_{\mathrm{0}}^L(I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right){\Psi }^{\prime }\mathrm{\delta }{\mathrm{\Phi }}^{\prime }\mathrm{d}x\hfill \\ +\mathrm{0.5}{\int }_{\mathrm{0}}^L\rho {\dot{\alpha }}^{\mathrm{2}}(I_z-I_y)\sin \left(\mathrm{2}\mathit{\gamma }\right)\Psi \mathrm{\delta }\mathrm{\Phi }\mathrm{d}x\hfill \\ +E{\int }_{\mathrm{0}}^L(I_y{\cos }^{\mathrm{2}}\mathit{\gamma }+I_z{\sin }^{\mathrm{2}}\mathit{\gamma }){\mathrm{\Phi }}^{\prime }\mathrm{\delta }{\mathrm{\Phi }}^{\prime }\mathrm{d}x\hfill \\ +{\kappa }_{z}G{\int }_{\mathrm{0}}^{L}A(W^{\prime }+\mathrm{\Phi })\mathrm{\delta }\mathrm{\Phi }\mathrm{d}x+k_{\mathrm{r}y}\mathrm{\Phi }\mathrm{\delta }\mathrm{\Phi }{|}_{x=\mathrm{0}}\hfill \end{array}\right\}=\mathrm{0}

Partial differential equations of motion in the rotz^d-direction:

(13)\left\{\begin{array}{l}{\int }_{\mathrm{0}}^L\rho (I_z{\cos }^{\mathrm{2}}\mathit{\gamma }+I_y{\sin }^{\mathrm{2}}\mathit{\gamma })\ddot{\Psi }\mathrm{\delta }\Psi \mathrm{d}x\hfill \\ -{\int }_{\mathrm{0}}^L\rho {\dot{\alpha }}^{\mathrm{2}}(I_z{\cos }^{\mathrm{2}}\mathit{\gamma }+I_y{\sin }^{\mathrm{2}}\mathit{\gamma })\Psi \mathrm{\delta }\Psi \mathrm{d}x\hfill \\ +\mathrm{2}{\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\sin \mathit{\gamma }\cos \beta -I_z\cos \mathit{\gamma }\sin \beta )\dot{\mathrm{\Theta }}\mathrm{\delta }\Psi \mathrm{d}x\hfill \\ +{\int }_{\mathrm{0}}^L\rho \ddot{\alpha }(I_y\sin \mathit{\gamma }\cos \beta -I_z\cos \mathit{\gamma }\sin \beta )\mathrm{\Theta }\mathrm{\delta }\Psi \mathrm{d}x\hfill \\ +{\mathrm{0.5}\int }_{\mathrm{0}}^L\rho (I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right)\ddot{\mathrm{\Phi }}\mathrm{\delta }\Psi \mathrm{d}x\hfill \\ +{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}E\int }_{\mathrm{0}}^L(I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right){\mathrm{\Phi }}^{\prime }\mathrm{\delta }{\Psi }^{\prime }\mathrm{d}x\hfill \\ -\mathrm{0.5}{\int }_{\mathrm{0}}^L\rho {\dot{\alpha }}^{\mathrm{2}}(I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right)\mathrm{\Phi }\mathrm{\delta }\Psi \mathrm{d}x\hfill \\ +{\int }_{\mathrm{0}}^L\rho \ddot{\alpha }(I_y\sin \mathit{\gamma }\sin \beta +I_z\cos \mathit{\gamma }\cos \beta )\mathrm{\delta }\Psi \mathrm{d}x\hfill \\ +E{\int }_{\mathrm{0}}^L(I_z{\cos }^{\mathrm{2}}\mathit{\gamma }+I_y{\sin }^{\mathrm{2}}\mathit{\gamma }){\Psi }^{\prime }\mathrm{\delta }{\Psi }^{\prime }\mathrm{d}x\hfill \\ -{\kappa }_{y}G{\int }_{\mathrm{0}}^{L}A(V^{\prime }+\Psi )\mathrm{\delta }\Psi \mathrm{d}x+k_{\mathrm{r}z}\Psi \mathrm{\delta }\Psi {|}_{x=\mathrm{0}}\hfill \end{array}\right\}=\mathrm{0}

The assumed mode method is utilized to discretize Eqs. (8)–(13) with the first N modes. By introducing canonical coordinates U_i(t), V_i(t), W_i(t), Q_i(t), F_i(t), and Y_i(t), the displacements U, V, W, Q, F, and Y can be expressed as follows:

(14)$$\begin{gathered} U={\displaystyle \sum _{i=\mathrm{1}}^N{\phi }_{\mathrm{1}i}\left(x\right)}{U}_i\left(t\right),V={\displaystyle \sum _{i=\mathrm{1}}^N{\phi }_{\mathrm{2}i}\left(x\right)}{V}_i\left(t\right), \\ W={\displaystyle \sum _{i=\mathrm{1}}^N{\phi }_{\mathrm{3}i}\left(x\right)}{W}_i\left(t\right),\mathrm{\Theta }={\displaystyle \sum _{i=\mathrm{1}}^N{\phi }_{\mathrm{4}i}\left(x\right)}{\mathrm{\Theta }}_i\left(t\right), \\ \mathrm{\Phi }={\displaystyle \sum _{i=\mathrm{1}}^N{\phi }_{\mathrm{5}i}\left(x\right)}{\mathrm{\Phi }}_i\left(t\right),\mathrm{\Psi }={\displaystyle \sum _{i=\mathrm{1}}^N{\mathit{\varphi }}_{\mathrm{6}i}\left(x\right)}{\Psi }_i\left(t\right), \end{gathered}$$

where Φ_ji(x) (j = 1, 2, …, 6; i = 1, 2, …, N) is the ith mode shape in the axial, flapwise, chordwise, and torsional rotation in the x- and y-axes and rotation in the z-axis. The corresponding expressions are as follows:

(15)\{\begin{array}{l}{\phi }_{\mathrm{1}i}=\cos \left[{\beta }_{\mathrm{a}i}\left(L-x\right)\right],\hfill \\ {\phi }_{\mathrm{2}i}=\cosh \left[{\beta }_{\mathrm{f}i}\left(L-x\right)\right]+\cos \left[{\beta }_{\mathrm{f}i}\left(L-x\right)\right]+{\xi }_{\mathrm{f}i}\left\{\sinh \left[{\beta }_{\mathrm{f}i}\left(L-x\right)\right]+\sin \left[{\beta }_{\mathrm{f}i}\left(L-x\right)\right]\right\},\hfill \\ {\phi }_{\mathrm{3}i}=\cosh \left[{\beta }_{\mathrm{c}i}\left(L-x\right)\right]+\cos \left[{\beta }_{\mathrm{c}i}\left(L-x\right)\right]+{\xi }_{\mathrm{c}i}\left\{\sinh \left[{\beta }_{\mathrm{c}i}\left(L-x\right)\right]+\sin \left[{\beta }_{\mathrm{c}i}\left(L-x\right)\right]\right\},\hfill \\ {\phi }_{\mathrm{4}i}=\cos \left[{\beta }_{\mathrm{t}i}\left(L-x\right)\right],\hfill \\ {\phi }_{\mathrm{5}i}={\displaystyle \frac{{{\phi }^{\prime }}_{\mathrm{3}i}}{i\mathrm{\pi }}},\hfill \\ {\phi }_{\mathrm{6}i}={\displaystyle \frac{{{\phi }^{\prime }}_{\mathrm{2}i}}{i\mathrm{\pi }}},\hfill \end{array}

where β_ai, β_fi, β_ci, and β_ti are the ith (i = 1, 2, …, N) roots of the characteristic equations in Eqs. (16a)–(16d).

(16a){\beta }_{\mathrm{a}}L\tan \left({\beta }_{\mathrm{a}}L\right)={\displaystyle \frac{k_{x}L}{E{A|}_{x=\mathrm{0}}}},

(16b)\left\{\begin{array}{l}{E}^{\mathrm{2}}{I}_z^{\mathrm{2}}{\left({\beta }_{\mathrm{f}}L\right)}^{\mathrm{4}}\left[\mathrm{1}-\cosh \left({\beta }_{\mathrm{f}}L\right)\cos \left({\beta }_{\mathrm{f}}L\right)\right]\hfill \\ -E{I}_z{\left({\beta }_{\mathrm{f}}L\right)}^{\mathrm{3}}L{k}_{\mathrm{r}z}\left[\cosh \left({\beta }_{\mathrm{f}}L\right)\sin \left({\beta }_{\mathrm{f}}L\right)+\sinh \left({\beta }_{\mathrm{f}}L\right)\cos \left({\beta }_{\mathrm{f}}L\right)\right]\hfill \\ +E{I}_z\left({\beta }_{\mathrm{f}}L\right)L^{\mathrm{3}}{k}_y\left[\sinh \left({\beta }_{\mathrm{f}}L\right)\cos \left({\beta }_{\mathrm{f}}L\right)-\cosh \left({\beta }_{\mathrm{f}}L\right)\sin \left({\beta }_{\mathrm{f}}L\right)\right]\hfill \\ +k_y{k}_{\mathrm{r}z}{L}^{\mathrm{4}}\left[\mathrm{1}+\cosh \left({\beta }_{\mathrm{f}}L\right)\cos \left({\beta }_{\mathrm{f}}L\right)\right]\hfill \end{array}\right\}=0，

(16c)\left\{\begin{array}{l}{E}^{\mathrm{2}}{I}_y^{\mathrm{2}}{\left({\beta }_{\mathrm{c}}L\right)}^{\mathrm{4}}\left[\mathrm{1}-\cosh \left({\beta }_{\mathrm{c}}L\right)\cos \left({\beta }_{\mathrm{c}}L\right)\right]\hfill \\ -E{I}_z{\left({\beta }_{\mathrm{c}}L\right)}^{\mathrm{3}}L{k}_{\mathrm{r}y}\left[\cosh \left({\beta }_{\mathrm{c}}L\right)\sin \left({\beta }_{\mathrm{c}}L\right)+\sinh \left({\beta }_{\mathrm{c}}L\right)\cos \left({\beta }_{\mathrm{c}}L\right)\right]\hfill \\ +E{I}_y\left({\beta }_{\mathrm{c}}L\right)L^{\mathrm{3}}{k}_z\left[\sinh \left({\beta }_{\mathrm{c}}L\right)\cos \left({\beta }_{\mathrm{c}}L\right)-\cosh \left({\beta }_{\mathrm{c}}L\right)\sin \left({\beta }_{\mathrm{c}}L\right)\right]\hfill \\ +k_z{k}_{\mathrm{r}y}{L}^{\mathrm{4}}\left[\mathrm{1}+\cosh \left({\beta }_{\mathrm{c}}L\right)\cos \left({\beta }_{\mathrm{c}}L\right)\right]\hfill \end{array}\right\}=0，

(16d){\beta }_{\mathrm{t}}L\tan \left({\beta }_{\mathrm{t}}L\right)={\displaystyle \frac{k_{\mathrm{r}x}L}{GJ}}.

In addition, the coefficients x_fi and x_ci in Eq. (15) can be calculated using the two following formulas:

(17a){\xi }_{\mathrm{f}i}=-{\displaystyle \frac{E{I}_z{\beta }_{\mathrm{f}i}\left[\cosh \left({\beta }_{\mathrm{f}i}L\right)-\cos \left({\beta }_{\mathrm{f}i}L\right)\right]+k_{\mathrm{r}\mathit{z}}\left[\sinh \left({\beta }_{\mathrm{f}i}L\right)-\sin \left({\beta }_{\mathrm{f}i}L\right)\right]}{E{I}_z{\beta }_{\mathrm{f}i}\left[\sin \left({\beta }_{\mathrm{f}i}L\right)-\sin \left({\beta }_{\mathrm{f}i}L\right)\right]+k_{\mathrm{r}\mathit{z}}\left[\cosh \left({\beta }_{\mathrm{f}i}L\right)+\cos \left({\beta }_{\mathrm{f}i}L\right)\right]}},

(17b){\xi }_{\mathrm{c}i}=-{\displaystyle \frac{E{I}_y{\beta }_{\mathrm{c}i}\left[\cosh \left({\beta }_{\mathrm{c}i}L\right)-\cos \left({\beta }_{\mathrm{c}i}L\right)\right]+k_{\mathrm{r}y}\left[\sinh \left({\beta }_{\mathrm{c}i}L\right)-\sin \left({\beta }_{\mathrm{c}i}L\right)\right]}{E{I}_y{\beta }_{\mathrm{c}i}\left[\sin \left({\beta }_{\mathrm{c}i}L\right)-\sin \left({\beta }_{\mathrm{c}i}L\right)\right]+k_{\mathrm{r}y}\left[\cosh \left({\beta }_{\mathrm{c}i}L\right)+\cos \left({\beta }_{\mathrm{c}i}L\right)\right]}}.

Substituting Eqs. (14), (15), and (17) into Eqs. (8)–(13), the equations of motion of the model can be expressed as follows:

(18)\mathit{M}\ddot{\mathit{q}}+\mathit{C}\dot{\mathit{q}}+\left({\mathit{K}}_{\mathrm{e}}+{\mathit{K}}_{\mathrm{c}}+{\mathit{K}}_{\mathrm{s}}+{\mathit{K}}_{\mathrm{acc}}\right)\mathit{q}=\mathit{F},

where M, C, K_e, K_c, K_s, and K_acc are mass, Coriolis force, structural stiffness, centrifugal stiffening, spin softening, and angular acceleration-induced stiffness matrices, respectively, and q and F are canonical coordinates and external force vector, respectively.

Model verification

In this section, the supporting stiffness of the proposed model (i.e., k_x, k_y, k_z, k_rx, k_ry, k_rz) is first calibrated via the experimental tests in combination with the genetic algorithm embedded in MATLAB toolbox. Then, the dynamic frequency characteristics of the calibrated semi-analytical model are verified based on the results obtained from the ANSYS models. Here, the specimen is made out of 2A12-T4 aluminum alloy. The corresponding material and geometrical parameter settings are listed as below:

Material properties: Young’s modulus E = 71.7 GPa; Density r = 2770 kg/m³; and Poisson’s ratio ν= 0.33.

Geometrical properties: Beam length L = 160 mm; Beam root width b₀= 25.2 mm; Beam root thickness h₀= 4.2 mm; Disk radius R_d = 150 mm; Tapered ratio of thickness t_h = 0.4048; Tapered ratio of width t_b = 0.3056; and Shear factors k_y = k_z = 5/6.

Calibration experiment on supporting stiffness

Before calibrating the elastic supporting stiffness, hammering and sweep-frequency tests are first conducted, as shown in Fig. 2. In the hammering test, additional testing instruments such as a force-hammer (PCB086C01) and a laser vibrometer (Polytec PDV-10) are adopted (see Fig. 2(a)). In the sweep-frequency test, a piezoelectric accelerometer (CA-YD-125, 1.5 g) and a shaker table (EDM-5000, 5–2500 Hz) are utilized (see Fig. 2(b)). Figure 2(c) shows the LMS front controller of a data acquisition system.

Based on the test system shown in Fig. 2, the first three-order natural frequencies of the tested specimen under β(L) = 45° are 164.3, 578.5, and 1149 Hz for the hammering test (see Fig. 3(a)) and 165, 580.5, and 1150 Hz for the frequency-sweep test (see Fig. 3(b)). Specifically, the three sections of frequency sweep are adopted in the sweep-frequency test for the sake of reducing the frequency-sweep time and exciting the first three-order natural frequencies as far as possible (see Figs. 2(b) and 3(b)). As a whole, the results obtained from both tests agree well with each other. In particular, the first three-order natural frequencies obtained from the hammering test are utilized to establish the optimized objective function Z shown in Eq. (19).

(19)Z=\sqrt{{\left(f_{\mathrm{n1}}\left(\mathit{k}\right)-\mathrm{164.3}\right)}^{\mathrm{2}}+{\left(f_{\mathrm{n2}}\left(\mathit{k}\right)-\mathrm{578.5}\right)}^{\mathrm{2}}+{\left(f_{\mathrm{n3}}\left(\mathit{k}\right)-\mathrm{1149}\right)}^{\mathrm{2}}},

where f_ni (i= 1, 2, 3) is the state variables, and k = [k_x, k_y, k_z, k_rx, k_ry, k_rz] is a vector related to the six design variables. Considering that the beam is clamped free, the supporting stiffness should be not too small. On this basis, the variation range of the supporting stiffness is all presupposed as [1 × 10⁶, 2 × 10⁸], and then the genetic algorithm is adopted to find the optimal stiffness combinations minimizing Z. The flow chart of the genetic algorithm adopted is shown in Fig. 4. In this work, the optimal k_x, k_y, k_z, k_rx, k_ry, and k_rz are finally determined as follows: k_x = 7.26 × 10⁶ N/m, k_y = 1.4 × 10⁶ N/m, k_z = 2 × 10⁷ N/m, k_rx = 3.7 × 10⁶ N·m/rad, k_ry = 3.8 × 10⁶ N·m/rad, and k_rz = 2.5 × 10⁶ N·m/rad. The first six characteristic roots (see Eqs. (16a)–(17b)) related to the optimal supporting stiffness are listed in as below:

Fig.2 Test rig: (a) Hammer test, (b) sweep test, and (c) data acquisition system.

Full size|PPT slide

Fig.3 Experimental results under γ(L) = 45°: (a) Frequency response function (FRF) of velocity obtained from hammering test and (b) spectrum cascades obtained from the sweep-frequency test. SFR: Sweep-frequency range; 1, 2, and 3: The first, the second, and the third sections of frequency sweep, respectively.

Full size|PPT slide

Fig.4 Flow chart of the genetic algorithm.

Full size|PPT slide

Support stiffness:

k_x = 7.26 ×10⁶ N/m;

k_y = 1.4 ×10⁶ N/m;

k_z = 2 ×10⁷ N/m;

k_rx= 3.7×10⁶ N·m/rad;

k_ry = 3.8 ×10⁶ N·m/rad;

k_rz = 2.5 ×10⁶ N·m/rad.

Characteristic root:

β_a=[0.3815, 3.1895, 6.3074, 9.4410, 12.5785, 15.7177];

β_f=[1.8681, 4.4769, 6.8342, 9.1140, 11.9561, 15.0046];

β_c=[1.8565, 4.1662, 6.1601, 8.8109, 11.8429, 14.9494];

β_t=[1.5708, 4.7123, 7.8538, 10.9953, 14.1368, 17.2783];

ξ_c=[−0.7226, −1.0427, −0.9963, −1.001, −1, −1];

ξ_f=[−0.7295, −1.0273, −0.9970, −1.001, −1, −1].

In Table 1, the results obtained from the proposed model with elastic and fixed supports are compared with those obtained from the experimental tests. As shown in Table 1, the results, especially for the higher-order modes under elastic-support boundary, are closer to the experimental results than those under fixed support boundary, which also indicates the impracticability of applying the infinite rigidity to the cantilevered end of the pre-twisted beam based on the fixture shown in Figs. 2(a) and 2(b).

Tab.1 First three-order natural frequency obtained from the proposed model and experiment

Mode	Frequencies obtained from hammering test (benchmark)/Hz	Frequencies obtained from sweep-frequency test/Hz	Frequencies obtained from proposed model/Hz
			Fixed support	Elastic support
fn1	164.3	165.0 (0.426%)	160.4 (−2.374%)	159.6 (−2.861%)
fn2	578.5	580.5 (0.346%)	598.2 (3.405%)	580.9 (0.415%)
fn3	1149.0	1150.0 (0.087%)	1208.6 (5.187%)	1149.8 (0.070%)

Note: The values in the brackets are errors relative to the results obtained from the hammering test.

ANSYS verification

In Subsection 2.3.1, the natural characteristics obtained from the proposed model under rotating speed n = 0 r/min are calibrated based on the experiments. Here, the corresponding dynamic characteristics are verified against the results obtained in the ANSYS models. In this work, Beam188, Shell181, Matrix27, and Mass21 elements in ANSYS software are adopted to build the finite element model of the pre-twisted beam with elastic supports, as shown in Fig. 5. Here, the elastic supports are simulated via the Matrix27 element, and the corresponding matrix formulation is asymmetric. All the displacement components Q = [U, V, W, Q, F, Y]^T of the semi-analytical model are determined in o^dx^dy^dz^d, and the supporting stiffness k_x, k_y, k_z, k_rx, k_ry, and k_rz in the semi-analytical model are also defined in o^dx^dy^dz^d, whereas the finite element model is established in O^rX^rY^rZ^r. Therefore, to make the semi-analytical model and the semi-analytical model equivalent, some necessary stiffness matrix transformation from o^dx^dy^dz^d to OX^rY^rZ^r is needed and presented as follows:

(20)$$\begin{gathered} {\mathit{K}}_{\mathrm{Sup,}{O}^{\mathrm{r}}{X}^{\mathrm{r}}{Y}^r{Z}^{\mathrm{r}}}={\mathit{T}}_{\mathrm{1}}^{\mathrm{T}}\left[\begin{array}{cccccc}\hfill k_x\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill k_y\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill k_z\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill k_{\mathrm{r}x}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill k_{\mathrm{r}y}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill k_{\mathrm{r}z}\hfill \end{array}\right]{\mathit{T}}_{\mathrm{1}}, \\ {\mathit{T}}_{\mathrm{1}}=\left[\begin{array}{cccccc}\hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \cos {\beta }_{\mathrm{0}}\hfill & \hfill \sin {\beta }_{\mathrm{0}}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill -\sin {\beta }_{\mathrm{0}}\hfill & \hfill \cos {\beta }_{\mathrm{0}}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \cos {\beta }_{\mathrm{0}}\hfill & \hfill \sin {\beta }_{\mathrm{0}}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill -\sin {\beta }_{\mathrm{0}}\hfill & \hfill \cos {\beta }_{\mathrm{0}}\hfill \end{array}\right]·\left[-{\mathit{I}}_{\mathrm{6}\times \mathrm{6}},{\mathit{I}}_{\mathrm{6}\times \mathrm{6}}\right], \end{gathered}$$

where {\mathit{K}}_{\mathrm{Sup,}{O}^{\mathrm{r}}{X}^{\mathrm{r}}{Y}^{\mathrm{r}}{Z}^{\mathrm{r}}}^{{}_{}} is the support-stiffness matrix in O^rX^rY^rZ^r, T₁ is the transfer matrix from o^dx^dy^dz^d to O^rX^rY^rZ^r, and I_6×6 is the 6 × 6 identity matrix.

Fig.5 Finite element model of a pre-twisted beam with elastic-support boundary: (a) Beam188 and (b) Shell181.

Full size|PPT slide

In terms of the Shell181 model shown in Fig. 5(b), the left end should be rigid with the help of the “CERIG” command and Mass21 element in ANSYS to make the Matrix27 element available in simulating the elastic supports.

Figure 6 shows the first-five order dynamic frequency variations obtained from the proposed model, the Beam188 model, and the Shell181 model. Generally speaking, the results obtained from the three models agree well with each other. However, the torsional dynamic frequencies (see f_n5 in Fig. 6) obtained from the proposed model are closer to those obtained from the Shell181 model rather than the Beam188 model. This is because the effect of the pre-twisted angle on the torsional stiffness in the Beam188 model is ignored [⁴³]. In addition, the first-five order modal shapes obtained from the proposed model under n = 0 r/min are also plotted in Fig. 6.

Fig.6 First five-order dynamic frequencies and modal shapes.

Full size|PPT slide

Analysis of modal coupling among y^d, z^d, and rotx^d vibrations

To quantify the degree of coupling among y^d, z^d, and rotx^d vibration modes varying with β₀, γ(L), and n, the modal strain energy method is then adopted to determine the corresponding MSERs. The MSER in the x^d direction, that is, the axial direction, is inessential and ignored, which is because the x^d vibration mode is hard to be excited in most cases. In addition, the effect of \ddot{\alpha } in Eqs. (8)–(13) is also ignored in this work. The QR method is used to solve Eq. (18) in order to obtain the rth corresponding right eigenvector q_r with a MATLAB program, that is,

(21)\mathit{M}{\ddot{\mathit{q}}}_r+\mathit{C}{\dot{\mathit{q}}}_r+\left({\mathit{K}}_{\mathrm{e}}+{\mathit{K}}_{\mathrm{c}}+{\mathit{K}}_{\mathrm{s}}+{\mathit{K}}_{\mathrm{acc}}\right){\mathit{q}}_r=\mathrm{0},

where {\mathit{q}}_r=[V_{\mathrm{1}},V_{\mathrm{2}},\mathrm{...},V_N,W_{\mathrm{1}},W_{\mathrm{2}},\mathrm{...},W_N,{\mathrm{\Theta }}_{\mathrm{1}},{\mathrm{\Theta }}_{\mathrm{2}},\mathrm{...},{\mathrm{\Theta }}_N,{\mathrm{\Phi }}_{\mathrm{1}}, {{\mathrm{\Phi }}_{\mathrm{2}},\mathrm{...},{\mathrm{\Phi }}_N,{\Psi }_{\mathrm{1}},{\Psi }_{\mathrm{2}},\mathrm{...},{\Psi }_N]}^{\mathrm{T}}.

Then, the MSER in the c direction (c = y^d, z^d, or rotx^d) can be calculated using the following formula:

(22)\mathrm{MSER}={\displaystyle \frac{\Vert {\left({\mathit{q}}_r^{\chi }\right)}^{\mathrm{T}}{\left({\mathit{K}}_{\mathrm{e}}+{\mathit{K}}_{\mathrm{c}}+{\mathit{K}}_{\mathrm{s}}\right)}_{\chi }{\mathit{q}}_r^{\chi }\Vert }{{\displaystyle \sum _{\chi =y^{\mathrm{b}},z^{\mathrm{b}},\mathrm{rot}{x}^{\mathrm{b}}}^{}\Vert {\left({\mathit{q}}_r^{\chi }\right)}^{\mathrm{T}}{({\mathit{K}}_{\mathrm{e}}+{\mathit{K}}_{\mathrm{c}}+{\mathit{K}}_{\mathrm{s}})}_{\chi }{\mathit{q}}_r^{\chi }\Vert }}},

where {\mathit{q}}_r^{y^{\mathrm{d}}}={\left[V_{\mathrm{1}},V_{\mathrm{2}},\mathrm{...},V_N,{\Psi }_{\mathrm{1}},{\Psi }_{\mathrm{2}},\mathrm{...},{\Psi }_N\right]}^{\mathrm{T}}, {\mathit{q}}_r^{z^{\mathrm{d}}}=[W_{\mathrm{1}},W_{\mathrm{2}},\mathrm{...},{W_N,{\mathrm{\Phi }}_{\mathrm{1}},{\mathrm{\Phi }}_{\mathrm{2}},\mathrm{...},{\mathrm{\Phi }}_N]}^{\mathrm{T}}, {\mathit{q}}_r^{\mathrm{rot}{x}^{\mathrm{d}}}={\left[{\mathrm{\Theta }}_{\mathrm{1}},{\mathrm{\Theta }}_{\mathrm{2}},\mathrm{...},{\mathrm{\Theta }}_N\right]}^{\mathrm{T}}, and ({\mathit{K}}_{\mathrm{e}}+{{\mathit{K}}_{\mathrm{c}}+{\mathit{K}}_{\mathrm{s}})}_{\chi } represents the stiffness matrix in the c direction.

Effects of β₀ on MSERs

Under the premise of γ(L) = 0°, the MSERs in the y^d, z^d, and rotx^d directions under β₀ = 0°, 15°, 30°, and 45° are plotted in Fig. 7. As can be seen, the degree of coupling in the y^d, z^d, and rotx^d directions increase with increasing n and β₀ other than β₀ = 0°. The reason for this is due to the existence of Coriolis coupling, stiffness coupling, and mass coupling terms in Eqs. (9)–(13). Some main conclusions are made as follows:

(1) When β₀ = 0° and n ≠ 0 r/min, the coupling terms between y^d and z^d and rotx^d directions do not exist (see +\mathrm{0.5}{\displaystyle {\int }_{\mathrm{0}}^L\rho A{\dot{\alpha }}^{\mathrm{2}}\sin \left(\mathrm{2}{\beta }_{\mathrm{0}}\right)W\mathrm{\delta }V\mathrm{d}x} in Eq. (9); +\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\sin \mathit{\gamma }}{\displaystyle \cos \beta -I_z\cos \mathit{\gamma }\sin \beta )\dot{\mathrm{\Theta }}\mathrm{\delta }\Psi \mathrm{d}x},+\mathrm{0.5}{\displaystyle {\int }_{\mathrm{0}}^L\rho (I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right)\ddot{\mathrm{\Phi }}\mathrm{\delta }\Psi }\mathrm{d}x, +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle {E\int }_{\mathrm{0}}^L(I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right){\mathrm{\Phi }}^{\prime }\mathrm{\delta }{\Psi }^{\prime }\mathrm{d}x}, and -\mathrm{0.5}{\displaystyle {\int }_{\mathrm{0}}^L\rho {\dot{\alpha }}^{\mathrm{2}}\left(I_y-I_z\right)\sin \left(\mathrm{2}\mathit{\gamma }\right)}\mathrm{\Phi }\mathrm{\delta }\Psi \mathrm{d}x in Eq. (13); and Figs. 7(a), 7(b), and 7(d)), but the coupling term between z^d and rotx^d directions still exists, leading to the appearance of MSERs in both directions (see the Coriolis coupling term -\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \mathit{\gamma }}{\displaystyle ·\cos \beta +I_z\sin \mathit{\gamma }\sin \beta )\dot{\mathrm{\Phi }}\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x} in Eq. (11) and Figs. 7(c) and 7(e)). In particular, for the y^d-dominated modes, the MSERs in the z^d and rotx^d directions are zero due to the nonexistence of the displacement component in both directions, and vice versa (see Figs. 7(a), 7(b), and 7(d)). Similar patterns can also be found in the z^d-/rot^d-dominated modes (see Figs. 7(c) and 7(e)). Here, the modal shapes corresponding to Figs. 7(a)–7(e) under 0 r/min are plotted in Fig. 8.

(2) When β₀≠ 0° and n ≠ 0 r/min, the degree of coupling among y^d, z^d, and rotx^d directions is dependent on the Coriolis coupling term -\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \beta \sin \mathit{\gamma }-I_z\cos \mathit{\gamma }}{\displaystyle ·\sin \beta )\dot{\Psi }\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x} in Eq. (11) and the stiffness coupling term +\mathrm{0.5}{\displaystyle {\int }_{\mathrm{0}}^L\rho A{\dot{\alpha }}^{\mathrm{2}}\sin \left(\mathrm{2}{\beta }_{\mathrm{0}}\right)W\mathrm{\delta }V\mathrm{d}x} in Eq. (9). Obviously, the absolute values of these coupling terms increase with the increase of β₀ (β₀∈[0°, 45°]). In other words, the degree of coupling among the y^d, z^d, and rotx^d directions is also gradually strengthened. More specifically, in terms of the flapwise-dominated modes (see Figs. 8(a), 8(b), and 8(d)), the MSERs in z^d and rotx^d directions are maximum, but the MSER in the y^d direction is minimum when β₀ = 45°. Similar rules are also suitable for the analysis of the flapwise- and torsional-dominated modes (see Figs. 7(c), 7(e), 8(c), and 8(e)).

Fig.7 Effects of β₀ on y^d, z^d, and rotx^d strain energy ratios varying with n: (a) f_n1, (b) f_n2, (c) f_n3, (d) f_n4, and (e) f_n5. Note: The left column represents MSER in the y^d direction; the middle column represents MSER in the z^d direction; and the right column represents MSER in the rotx^d direction.

Full size|PPT slide

(3) A larger n leads to a stronger degree of coupling among y^d, z^d, and rotx^d directions for the same β₀ (Figs. 7(a)–7(e)). This can be easily illustrated via the aforementioned conclusions (1) and (2), and no more expatiation is presented here.

Fig.8 First five-order modal shapes under n = 0 r/min: (a) f_n1, (b) f_n2, (c) f_n3, (d) f_n4, and (e) f_n5.

Full size|PPT slide

Effects of γ(L) on MSERs

Under the premise of β₀ = 0°, the MSERs in the y^d, z^d, and rotx^d directions under γ(L) = 0°, 15°, 30°, and 45° are shown in Fig. 9. Here, the effect of γ(L) = 0° on MSERs has been explicitly analyzed in Subsection 3.1, and no more clarifications are made here. Obviously, the introduction of γ(L) automatically leads to the structural coupling between y^d and z^d directions (see the mass coupling term +0.5 ·{\displaystyle {\int }_{\mathrm{0}}^L\rho (I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right)\ddot{\Psi }\mathrm{\delta }\mathrm{\Phi }\mathrm{d}x} and structural stiffness coupling term +\mathrm{0.5}E{\displaystyle {\int }_{\mathrm{0}}^L(I_y-I_z)\sin \left(\mathrm{2}\mathit{\gamma }\right){\Psi }^{\prime }\mathrm{\delta }{\mathrm{\Phi }}^{\prime }\mathrm{d}x} in Eq. (12)). In addition, the introduction of n further causes the coupling among the y^d, z^d, and rotx^d directions (see the Coriolis coupling term +\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \mathit{\gamma }\cos \beta +I_z\sin \mathit{\gamma }\sin \beta )\dot{\mathrm{\Phi }}\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x}-\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \beta \sin \mathit{\gamma }-I_z\cos \mathit{\gamma }\sin \beta )\dot{\Psi }\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x} in Eq. (11) and stiffness coupling term +\mathrm{0.5}{\displaystyle {\int }_{\mathrm{0}}^L\rho {\dot{\alpha }}^{\mathrm{2}}\left(I_z-I_y\right)\sin \left(\mathrm{2}\mathit{\gamma }\right)\Psi \mathrm{\delta \Phi d}x} in Eq. (12)). Some main conclusions are summarized as follows:

(1) In terms of flapwise-dominated modes (see Figs. 6, 9(a), 9(b), and 9(d)), MSER decreases in the y^d direction and increases in the z^d direction with increasing γ(L), which indicates that the degree of coupling is increasingly strengthened. This can be attributed to the enhanced coupling effects such as mass coupling in Eq. (11) and stiffness coupling in Eqs. (11) and (12). In addition, the MSER in the rotx^d direction comes from the two coupling contributions in both y^d and z^d directions (see the Coriolis coupling term +\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \mathit{\gamma }\cos \beta +I_z\sin \mathit{\gamma }\sin \beta )\dot{\mathrm{\Phi }}\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x}-\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \beta \sin \mathit{\gamma }-I_z\cos \mathit{\gamma }\sin \beta )\dot{\Psi }\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x} in Eq. (11)). The coupling contribution coming from the y^d direction increases, but those coming from the z^d direction decrease with increasing γ(L), causing the non-monotonic change of rotx^d MSER under the combined effect of y^d and z^d MSER (see Figs. 9(a), 9(b), and 9(d)). Similar rules also exist in the flapwise-dominated modes (see Figs. 6 and 9(c)) and need not be repeated here.

Fig.9 Effects of γ(L) on y^d, z^d, and rotx^d strain energy ratios varying with n: (a) f_n1, (b) f_n2, (c) f_n3, (d) f_n4, and (e) f_n5. Note: The left column represents MSER in the y^d direction; the middle column represents MSER in the z^d direction; and the right column represents MSER in the rotx^d direction.

Full size|PPT slide

(2) In terms of torsional-dominated mode (see Figs. 6 and 9(e)), the MSER in the rotx^d direction decreases with increasing γ(L), whereas those in the y^d and z^d directions increase and decrease with increasing γ(L), respectively, which is mainly due to the strengthened coupling effect in the y^d direction (see the Coriolis coupling term −2·{\displaystyle {\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \beta \sin \mathit{\gamma }-I_z\cos \mathit{\gamma }\sin \beta )\dot{\Psi }\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x} in Eq. (11)) and the degraded coupling effect in the z^d direction (see the Coriolis coupling term -\mathrm{2}{\displaystyle {\int }_{\mathrm{0}}^L\rho \dot{\alpha }(I_y\cos \mathit{\gamma }\sin \beta +I_z\sin \mathit{\gamma }\sin \beta )} {\displaystyle \dot{\mathrm{\Phi }}\mathrm{\delta }\mathrm{\Theta }\mathrm{d}x} in Eq. (11)).

(3) For the flapwise-/chordwise-dominated mode (see Figs. 6 and 9(a)–9(d)), the MSERs in the y^d and z^d directions are complementary in the changing trend and close to 1 with increasing n for the same γ(L), whereas those for the torsional-dominated mode (see Figs. 6 and 9(e)) are the same in the changing trend.

Conclusions

In this work, Hamilton’s principle in combination with assumed mode method is adopted to establish the coupled flapwise–chordwise–axial–torsional dynamic model of a rotating double-tapered, pre-twisted, and inclined Timoshenko beam with elastically restrained root. Then, the proposed model is verified by comparing natural characteristics obtained from the experiment and ANSYS models using Beam188 and Shell181 elements. Finally, the MSER is adopted to quantify the degree of coupling among the flapwise, chordwise, and torsional directions under the effects of setting angle, pre-twisted angle, and rotating speed. Some conclusions are made as follows:

(1) Both flapwise–torsional coupling and chordwise–torsional coupling are determined via the Coriolis coupling term. The degrees of flapwise–torsional coupling and chordwise-torsional coupling increase and decrease with increasing setting and pre-twisted angles, respectively. The changing trend of MSER in the torsional direction is dependent on the combined effect of the flapwise–torsional and chordwise–torsional couplings, especially for flapwise-dominated or chordwise-dominated mode.

(2) The pre-twisted angle has more effect on the degree of flapwise–chordwise coupling than the setting angle, thereby inducing stronger mode coupling. The introduction of pure setting angle only causes one stiffness coupling term related to the rotating speed, whereas the introduction of the pure pre-twisted angle not only leads to the stiffness coupling term related to the rotating speed but also the structural stiffness and mass couplings independent on the rotating speed.

Acknowledgements

This project was supported by the National Natural Science Foundation (Grant Nos. 11972112 and 11772089), the Fundamental Research Funds for the Central Universities (Grant Nos. N170308028, N170306004, N2003014, and N180708009), and Liaoning Revitalization Talents Program (Grant No. XLYC1807008).

References

1 Lin S M. Dynamic analysis of rotating nonuniform Timoshenko beams with an elastically restrained root. Journal of Applied Mechanics, 1999, 66(3): 742–749

2 He Q, Xuan H J, Liu L L, et al. Perforation of aero-engine fan casing by a single rotating blade. Aerospace Science and Technology, 2013, 25(1): 234–241

3 Javdani S, Fabian M, Carlton J S, et al. Underwater free-vibration analysis of full-scale marine propeller using a fiber Bragg grating-based sensor system. IEEE Sensors Journal, 2016, 16(4): 946–953

4 Rezaei M M, Behzad M, Haddadpour H, et al. Development of a reduced order model for nonlinear analysis of the wind turbine blade dynamics. Renewable Energy, 2015, 76: 264–282

5 Rao J S, Carnegie W. Solution of the equations of motion of coupled-bending bending torsion vibrations of turbine blades by the method of Ritz–Galerkin. International Journal of Mechanical Sciences, 1970, 12(10): 875–882

6 Houbolt J C, Brooks G W. Differential equations of motion for combined flapwise bending, chordwise bending, and torsion of twisted nonuniform rotor blades. National Advisory Committee for Aeronautics, Technical Note 3905, 1957

7 Du H, Lim M K, Liew K M. A power series solution for vibration of a rotating Timoshenko beam. Journal of Sound and Vibration, 1994, 175(4): 505–523

8 Rao J S. Flexural vibration of pretwisted tapered cantilever blades. Journal of Engineering for Industry, 1972, 94(1): 343–346

9 Şakar G, Sabuncu M. Dynamic stability of a rotating asymmetric cross-section blade subjected to an axial periodic force. International Journal of Mechanical Sciences, 2003, 45(9): 1467–1482

10 Banerjee J R. Free vibration of centrifugally stiffened uniform and tapered beams using the dynamic stiffness method. Journal of Sound and Vibration, 2000, 233(5): 857–875

11 Zeng J, Ma H, Yu K, Rubbing response comparisons between single blade and flexible ring using different rubbing force models. International Journal of Mechanical Sciences, 2019, 164: 105164

12 Zeng J, Zhao C G, Ma H, Rubbing dynamic characteristics of the blisk-casing system with elastic supports. Aerospace Science and Technology, 2019, 95: 105481

13 Hodges D H, Dowell E H. Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades. NASA Technical Report NASA-TN-D-7818, A-5711, 1974

14 Zhu T L. The vibrations of pre-twisted rotating Timoshenko beams by the Rayleigh–Ritz method. Computational Mechanics, 2011, 47(4): 395–408

15 Ma H, Xie F T, Nai H Q, Vibration characteristics analysis of rotating shrouded blades with impacts. Journal of Sound and Vibration, 2016, 378: 92–108

16 Xie F T, Ma H, Cui C, Vibration response comparison of twisted shrouded blades using different impact models. Journal of Sound and Vibration, 2017, 397: 171–191

17 Sinha S K. Combined torsional-bending-axial dynamics of a twisted rotating cantilever Timoshenko beam with contact-impact loads at the free end. Journal of Applied Mechanics, 2007, 74(3): 505–522

18 Yang X D, Wang S W, Zhang W, et al. Dynamic analysis of a rotating tapered cantilever Timoshenko beam based on the power series method. Applied Mathematics and Mechanics, 2017, 38(10): 1425–1438

19 Şakar G, Sabuncu M. Buckling and dynamic stability of a rotating pretwisted asymmetric cross-section blade subjected to an axial periodic force. Finite Elements in Analysis and Design, 2004, 40(11): 1399–1415

20 Sabuncu M, Evran K. Dynamic stability of a rotating pre-twisted asymmetric cross-section Timoshenko beam subjected to an axial periodic force. International Journal of Mechanical Sciences, 2006, 48(6): 579–590

21 Subrahmanyam K B, Kaza K R V, Brown G V, et al. Nonlinear bending-torsional vibration and stability of rotating, pretwisted, preconed blades including Coriolis effects. In: Proceedings of Workshop on Dynamics and Aeroelastic Stability Modeling of Rotor Systems. Atlanta: NASA, 1986, NASA-TM-87207

22 Sina S A, Haddadpour H. Axial–torsional vibrations of rotating pretwisted thin walled composite beams. International Journal of Mechanical Sciences, 2014, 80: 93–101

23 Adair D, Jaeger M. Vibration analysis of a uniform pre-twisted rotating Euler–Bernoulli beam using the modified Adomian decomposition method. Mathematics and Mechanics of Solids, 2018, 23(9): 1345–1363

24 Oh Y, Yoo H H. Vibration analysis of a rotating pre-twisted blade considering the coupling effects of stretching, bending, and torsion. Journal of Sound and Vibration, 2018, 431: 20–39

25 Lee J W, Lee J Y. Development of a transfer matrix method to obtain exact solutions for the dynamic characteristics of a twisted uniform beam. International Journal of Mechanical Sciences, 2016, 105: 215–226

26 Banerjee J R. Development of an exact dynamic stiffness matrix for free vibration analysis of a twisted Timoshenko beam. Journal of Sound and Vibration, 2004, 270(1–2): 379–401

27 Subrahmanyam K B, Kaza K R V. Vibration and buckling of rotating, pretwisted, preconed beams including Coriolis effects. Journal of Vibration and Acoustics, 1986, 108(2): 140–149

28 Hashemi S M, Richard M J. Natural frequencies of rotating uniform beams with Coriolis effects. Journal of Vibration and Acoustics, 2001, 123(4): 444–455

29 Banerjee J R, Kennedy D. Dynamic stiffness method for inplane free vibration of rotating beams including Coriolis effects. Journal of Sound and Vibration, 2014, 333(26): 7299–7312

30 Oh S Y, Song O, Librescu L. Effects of pretwist and presetting on coupled bending vibrations of rotating thin-walled composite beams. International Journal of Solids and Structures, 2003, 40(5): 1203–1224

31 Latalski J, Warminski J, Rega G. Bending–twisting vibrations of a rotating hub–thin-walled composite beam system. Mathematics and Mechanics of Solids, 2017, 22(6): 1303–1325

32 Ondra V, Titurus B. Free vibration analysis of a rotating pre-twisted beam subjected to tendon-induced axial loading. Journal of Sound and Vibration, 2019, 461: 114912

33 Zeng J, Chen K K, Ma H, Vibration response analysis of a cracked rotating compressor blade during run-up process. Mechanical Systems and Signal Processing, 2019, 118: 568–583

34 Sun Q, Ma H, Zhu Y P, Comparison of rubbing induced vibration responses using varying-thickness-twisted shell and solid-element blade models. Mechanical Systems and Signal Processing, 2018, 108: 1–20

35 Ma H, Wang D, Tai X Y, Vibration response analysis of blade-disk dovetail structure under blade tip rubbing condition. Journal of Vibration and Control, 2017, 23(2): 252–271

36 Lin S M. The instability and vibration of rotating beams with arbitrary pretwist and an elastically restrained root. Journal of Applied Mechanics, 2001, 68(6): 844–853

37 Lin S M, Wu C T, Lee S Y. Analysis of rotating nonuniform pretwisted beams with an elastically restrained root and a tip mass. International Journal of Mechanical Sciences, 2003, 45(4): 741–755

38 Lee S Y, Lin S M, Wu C T. Free vibration of a rotating non-uniform beam with arbitrary pretwist, an elastically restrained root and a tip mass. Journal of Sound and Vibration, 2004, 273(3): 477–492

39 Choi S T, Chou Y T. Vibration analysis of elastically supported turbomachinery blades by the modified differential quadrature method. Journal of Sound and Vibration, 2001, 240(5): 937–953

40 Bambill D V, Rossit C A, Rossi R E, Transverse free vibration of non uniform rotating Timoshenko beams with elastically clamped boundary conditions. Meccanica, 2013, 48(6): 1289–1311

41 Digilov R M, Abramovich H. The impact of root flexibility on the fundamental frequency of a restrained cantilever beam. International Journal of Mechanical Engineering Education, 2017, 45(2): 184–193

42 Hodges D H. Torsion of pretwisted beams due to axial loading. Journal of Applied Mechanics, 1980, 47(2): 393–397

43 Zeng J, Ma H, Yu K, Coupled flapwise-chordwise-axial-torsional dynamic responses of rotating pre-twisted and inclined cantilever beams subject to the base excitation. Applied Mathematics and Mechanics, 2019, 40(8): 1053–1082