Nomenclature

A, cross-sectional area;

A₀, amplitude of the base excitation;

OXYZ, global coordinate system;

OX^rY^rZ^r, rotating coordinate system;

o^bx^by^bz^b, local coordinate system;

o^ex^ey^ez^e, element coordinate system;

A_i, coordinate-transformation matrix;

b, beam width at the arbitrary section;

b₀, beam width at the root section;

c, constant;

C, Rayleigh damping;

E, Young's modulus;

F_b, base excitation force in the x-direction;

F_c, centrifugal force;

F^e, element nodal force vector in o^ex^ey^ez^e;

F_r^e, element nodal force vector in OX^rY^rZ^r;

F, global nodal force vector;

f_e, excitation frequency;

G, shear modulus;

f_i, concerned end frequency;

G^e, element Coriolis matrix in o^ex^ey^ez^e;

G_r^e, element Coriolis matrix in OX^rY^rZ^r;

G, global Coriolis matrix;

h, beam thickness at the arbitrary section;

h₀, beam thickness at the root section;

I_x, area moment of the cross-section inertia about ox₃;

I_y, area moment of the cross-section inertia about oy₃;

I_z, area moment of the cross-section inertia about oz₃;

J, torsional moment of the inertia of a pre-twisted rectangular beam;

J_s, torsional moment of the inertia of a straight rectangular beam;

J_a, additional torsional moment of the inertia caused by γ(L);

K_e^e, element structural stiffness matrix in o^ex^ey^ez^e;

K_c^e, element centrifugal stiffening matrix in o^ex^ey^ez^e when =60r·min^-1;

K_s^e, spin softening matrix in o^ex^ey^ez^e when =60r·min^-1;

K_{e, r}^e, element structural stiffness matrix in OX^rY^rZ^r;

K_{c, r}^e, element centrifugal stiffening matrix in OX^rY^rZ^r when =60r·min^-1;

K_{s, r}^e, element spin softening matrix in OX^rY^rZ^r when =60r·min^-1;

K_e, global structural stiffness matrix;

K_c, global centrifugal stiffening matrix;

K_s, global spin softening matrix;

K_acc^e, stiffness matrix related to in o^ex^ey^ez^e;

K_{acc, r}^e, stiffness matrix related to in OX^rY^rZ^r;

K_acc, global stiffness matrix related to ;

L, beam length;

l_k, kth element length;

M^e, element mass matrix in o^ex^ey^ez^e;

M_r^e, element mass matrix in OX^rY^rZ^r;

M, global mass matrix;

N, total number of the beam elements;

n, rotation speed;

R_d, disk radius;

R_k, x coordinate of o^e in OX^rY^rZ^r;

r_P, coordinate vector of the arbitrary point Pin OXYZ;

, velocity vector of the arbitrary point Pin OXYZ;

T_k, kinetic energy of the kth rotating Timoshenko beam element;

T, transformation matrix from OX^rY^rZ^r to o^ex^ey^ez^e;

t₁, starting moment of the variation;

t₂, ending moment of the variation;

U_k, potential energy of the kth rotating Timoshenko beam element;

u, v, w, linear displacements of the centroid of an arbitrary beam cross-section along the positive x^e-, y^e-, and z^e-axes;

, velocities of the centroid of an arbitrary beam cross-section along the positive x^e-, y^e-, and z^e-axes;

, accelerations of the centroid of an arbitrary beam cross-section along the positive x^e-, y^e- and z^e-axes;

u', v', w', derivatives of u, v, and w versus the coordinate x;

X_u^r, base excitation along the X^r-direction;

x_k, x coordinate of o^e in o^bx^by^bz^b;

(X_P, Y_P, Z_P), coordinate components of Point P along the X-, Y-, and Z-directions;

(x, y, z), coordinate components of Point Palong the x₃-, y₃-, and z₃-directions;

Y_I, first-order flapwise frequency;

Y_II, second-order flapwise frequency;

Z_I, first-order chordwise frequency;

Z_II, second-order chordwise frequency.

Greek symbols

α, rotating angle;

, angular velocity;

, angular acceleration;

α₁, β₁, Rayleigh damping coefficients;

β₀, initial setting angle;

β_k, setting angle of the kth Timoshenko beam element;

γ_k, pre-twisted angle of the kth Timoshenko beam element;

γ(L), pre-twisted angle at L;

δ, variational symbol;

δ^e, element nodal displacement vector in o^ex^ey^ez^e;

δ_r^e, element nodal displacement vector in OX^rY^rZ^r;

ζ, correction factor;

θ, angular displacement with respect to the x^e-axis;

, velocity with respect to the x^e-axis;

, acceleration with respect to the x^e-axis;

θ', derivative of θ versus the coordinate x;

θ_I, first-order torsional frequency;

, shear coefficient in the y₃-direction;

, shear coefficient in the z₃-direction;

ξ₁, ξ₂, modal damping ratios;

ρ, density;

τ_b, breadth taper;

τ_h, thickness taper;

υ, Poisson's ratio;

φ, angular displacement with respect to the z^e-axis;

, velocity with respect to the z^e-axis;

, acceleration with respect to the ze axis;

φ', derivative of φ versus the coordinate x;

χ, symbol representing global matrices;

χ_k^e, symbol representing element matrices;

ϕ, angular displacement with respect to the y^e-axis;

, velocity with respect to the y^e-axis;

, acceleration with respect to the y^e-axis;

ϕ', derivative of ϕ versus the coordinate x.

1 Introduction

Rotating beam-like structures are widely spread over many engineering applications such as compressor blades, helicopter blades, wind blades, and propellers^[1-2]. In general, almost all of the structural vibration problems such as fatigue failure, flutter, and stability are closely related to the natural frequencies of the structures^[3-5]. Therefore, it is essential to have a thorough understanding of their dynamic characteristics for engineering design purposes^[6].

The equations of motion of rotating beams have been studied extensively^[1-10]. Unlike non-rotating beams, rotating beams have additional Coriolis, centrifugal stiffening, and spin softening effects, which will result in the couplings of structural modes among different vibration directions^[9-10]. In Ref.[11], the coupled bending-bending equations of motion of a pre-twisted and tapered Euler beam with a fluctuant angular velocity were determined with the Hamilton principle, while the coupled bending-bending equations of motion of a rotating composite thin-walled beam subject to thermal stress field, humidity stress field, and aerodynamic forces were elaborately derived in Ref.[12]. The coupled flapwise-chordwise-axial equations of motion of the rotating uniform Timoshenko beam with a variable speed were established in Refs.[13] and[14]. In Refs.[5], [15], and [16], the effects of spin softening and Coriolis on the dynamics of a rotating cantilever beam were excluded when its coupled flapwise-chordwise-torsional equations of motion were derived. Such effects were also considered in the free vibration analyses of the rotating Timoshenko beam featuring flapwise-chordwise-torsional couplings in Ref.[17] and flapwise-axial-torsional couplings in Ref.[18]. In Refs.[2], [6], and [19], although the coupled flapwise-chordwise-axial-torsional equations of motion of a rotating Timoshenko beam were established, only the theoretical derivation was carried out. Some necessary verification works were deficient in Ref.[19], while the effects of centrifugal stiffening on the torsional behaviors of the beam were ignored in Refs.[2] and [6].

The equations of motion mentioned before^[1-19] are partial differential equations, and their closed-form solutions are not available in most cases^{[2, 20]}. To solve these partial differential equations, lots of numerical algorithms emerge, e.g., the dynamic stiffness method^[21-23], the transfer matrix method^[24], the differential transformation method^[25], the Galerkin method^{[13, 26]}, the Rayleigh-Ritz method^[18], the power series method^[27], the Chebyshev-Ritz method^[28], and the finite element (FE) method^[29-30]. Generally speaking, the FE method is more suitable for handling problems such as complicated geometrical configurations, loads, materials, and complicated boundary conditions than the other methods^[31-32].

The literature survey presented here indicates that little emphasis has been placed on the coupled flapwise-chordwise-axial-torsional dynamic modeling of rotating beams with the consideration of all rotating effects, i.e., the effects of Coriolis force, centrifugal stiffening, spin softening, and angular acceleration. Moreover, the modal characteristics of rotating beams are the concern of most scholars, and little attention has been paid to the coupled flapwise-chordwise-axial-torsional vibration responses under the base excitation. Aiming to these deficiencies, a rotating pre-twisted and inclined cantilever beam model (RPICBM) is established with the Hamilton principle and the FE method (see Section 2). Based on the proposed model shown in Section 2, the effects of the setting angle and the pre-twisted angle on the modal characteristics and the base excitation-induced vibration responses are analyzed theoretically (see Section 3). Finally, some conclusions are made in Section 4.

2 FE modeling of the RPICBM subject to the base excitation

An engineering example of a bladed disk structure subject to the base excitation is shown in Fig. 1(a), which consists of two first-stage fan blades, a rigid disk, and elastic supports^[33]. In Fig. 1(a), the intersection angle between the unreleased blade and the released blade is 180°. Once the test rig is rotating and reaches a certain rotation speed, one of the two first-stage fan blades will be released from the rigid disk under the effect of the centrifugal force, which leads to the unbalanced bladed disk system. Then, the unbalanced bladed disk system will be whirling. Meanwhile, the motion of the unreleased blade can be viewed as the resultant motion of the rigid disk whiling and blade rotating^[34], i.e., the dynamics of the unreleased rotating blade can be degenerated into analyzing the dynamic characteristics of a rotating cantilever blade subject to the base excitation^[35].

Based on the aforementioned engineering example, a similar mathematical model is simplified and adopted in this paper, i.e., an FE model of the RPICBM subject to the base excitation is pre-supposed and then established in this paper via combining the Hamilton principle with the FE method (see Fig. 1(c)). Figure 1(c) shows that the RPICBM is discretized into a series of inclined straight rotating Timoshenko beam segments, which are so-called Timoshenko beam elements. In this section, emphasis will be placed on the derivation of the partial differential equations and the element matrices of rotating Timoshenko beams. It is worth noting that the rotating effects such as Coriolis, centrifugal stiffening, spin softening, and angular acceleration are all included in the flapwise-, chordwise-, axial-, and torsional-directions. Finally, the proposed RPICBM is verified via comparisons of the natural frequencies and the base excitation-induced vibration responses among literature and FE models by using Beam188 and Shell181 elements in ANSYS.

2.1 Establishment of the partial differential equations

The kth Timoshenko beam element is taken as an example to establish the corresponding coupled flapwise-chordwise-axial-torsional equations of motion (see Fig. 1(c)). In Fig. 1(c), OXYZ, OX^rY^rZ^r, o^bx^by^bz^b, o^ex^ey^ez^eand are the global coordinate system, the rotating coordinate system, the local coordinate system, and the element coordinate system, respectively, and α is the rotation angle. In Fig. 1(b), three coordinates of the second-kind Euler angles θ, ϕ, and φ are adopted to describe the angular displacements of an arbitrary rigid beam cross-section in o^ex^ey^ez^e, β₀ is the initial setting angle, and γ_k is the pre-twisted angle of the kth Timoshenko beam element. Here,

is defined as the setting angle of the kth Timoshenko beam element. Let

where R_k is the x-coordinate of o^e in OX^rY^rZ^r, R_d is the disk radius, x_k is the x-coordinate of o^e in o^bx^by^bz^b, and X_u^r is the base excitation along the X^r-direction. Then, the coordinates of an arbitrary point P on the beam cross-section in OXYZ can be expressed as follows:

(1)

where (x, y, z) is the coordinates of Point P in ox₃y₃z₃ (see Fig. 1(b)). u, v, and w are the linear displacements of the centroid of an arbitrary beam cross-section along the positive x^e-, y^e-, and z^e-axes, respectively. A₁, A₂, A₃, A₄, and A₅ are the coordinate-transformation matrices from ox₃y₃z₃ to OXYZ, which can be expressed as follows:

(2)

In Eq.(2), the first-order approximation is adopted based on the small deformation assumption, i.e.,

According to Eqs.(1) and (2), the kinetic energy T_k of the kth rotating Timoshenko beam element can be calculated as follows:

(3)

where ρ, l_k, and A are the density, the kth element length, and the cross-sectional area, respectively. An over-dot in Eq.(3) represents the differentiation with respect to time.

According to Refs.[4], [17], [36], and [37], the potential energy U_k of the kth rotating Timoshenko beam element is shown as follows:

(4)

where E and G are the Young's modulus and the shear modulus, respectively. I_x, I_y, and I_z are the area moments of the inertia of the cross-section about ox₃, oy₃, and oz₃, respectively. J is the torsional moment of inertia, and its calculation is shown in Subsection 2.4. and are the shear coefficients in the y₃- and z₃-directions, respectively.

is the centrifugal force.

is the base excitation force in the x-direction.

is a constant. In addition, a prime in Eq.(4) represents the differentiation with respect to the coordinate x. In Eq.(4), it should be noted that the underlined term is the effect of the centrifugal stiffening on the torsional behavior, while the double underlined term is the contribution of the pre-twisted angle to the potential energy^[37].

Taking the functionals u, v, w, θ, ϕ, and φ as independent variables, the Hamilton principle is applied to the derived energy equation

and six partial differential equations are then determined. The x^e-axial displacement u is

(5)

The y^e-bending direction displacement v is

(6)

The z^e-bending direction displacement w is

(7)

The angular displacement θ with respect to the x^e-axis is

(8)

The angular displacement ϕ with respect to the y^e-axis is

(9)

The angular displacement φ with respect to the z^e-axis is

(10)

It should be noted that the terms including in Eqs.(5)–(10) are the angular acceleration-related terms.

2.2 Element matrices of a rotating Timoshenko beam

In order to solve Eqs.(5)–(10), the displacement field functions in the u-, v-, w-, θ-, ϕ-, and φ-directions are assumed, and the corresponding expressions are shown in Appendix A. Substituting Eqs.(A1)–(A4) into Eqs.(5)–(10), the following second-order ordinary differential equations are determined:

(11)

where M^e and k_e^e represent the mass and structural stiffness matrices, respectively. G^e, K_c^e, and K_s^e are the Coriolis, centrifugal stiffening, and spin softening matrices when the angular velocity =60 r·min^-1, respectively. is the stiffness matrix related to the angular acceleration . and indicate the nodal displacement and force vectors, respectively. More detailed information of these element matrices can refer to Appendix B.

2.3 Element matrix assembling

All the element matrices, i.e., M^e, G^e, k_e^e, K_c^e, K_s^e, K_acc^e, , and in Subsection 2.2 are determined in the element coordinate system o^ex^ey^ez^e. However, when these element matrices are assembled to establish the RPICBM, the corresponding element coordinate systems must be unified. Here, all the element matrices are assembled in the rotating coordinate system OX^rY^rZ^r, and the matrix transformation formulas are shown as follows:

(12)

where T is the transformation matrix from OX^rY^rZ^r to o^ex^ey^ez^e. and are the nodal displacement vectors in o^ex^ey^ez^e and OX^rY^rZ^r, respectively.

Based on Eqs.(11) and (12), the element matrices in OX^rY^rZ^r such as M_r^e, G_r^e, K_{e, r} ^e, K_{c, r}^e, K_{s, r}^e, K_{acc, r}^e, and F_r^e can be calculated as follows:

(13)

The global matrices of the RPICBM can be determined via . N is the total number of the beam elements.

2.4 Torsional stiffness of a pre-twisted rectangular beam

A detailed derivation process for the calculation of the torsional moment of inertia of a straight rectangular beam J_s with the consideration of the effect of the thickness-width ratio (h/b, see Fig.C1 in Appendix C) is shown in Appendix C. Here, the expressions for calculating J_s are directly shown as follows:

(14)

Furthermore, an additional torsional moment of inertia J_a caused by a pre-twisted angle γ(L) is determined as follows^[40]:

(15)

where γ(L) is the pre-twisted angle at the distance L (see Fig. 1).

In terms of Eqs.(14) and (15), J can be determined via J_s+J_a. In Fig. 2, the torsional stiffness GJ and the first-order torsional frequency θ_I obtained from this paper agree well with those obtained from Ref.[40], and this verifies the derived J.

2.5 Model validation

Based on Subsections 2.1–2.4, the FE model of the RPICBM is established in Fig. 3. In addition, some measures should be taken to validate the developed model. In terms of the effects of the profile tapper and the pre-twisted angle on the dynamic characteristics of a beam-like structure, various beam models have been proposed^{[11, 19, 21, 23, 30, 41]}. Here, a pre-twisted cantilever beam model in Ref.[41] is taken as an example to illustrate the availability of the developed model. The corresponding material and geometrical parameters are defined as follows^[41]:

where v is Poisson's ratio, L is the beam length, b₀ is the beam root width, h₀ is the beam root thickness, and τ_h is the thickness taper. In order to verify the proposed model, comparisons of the natural frequencies and the vibration responses are made among the proposed model, Ref.[41], and FE models by using Beam188 and Shell181 elements (see Figs. 4 and 5).

In Fig. 4, Y_I, Z_I, Y_II, and Z_II obtained from the proposed model, the FE models by using shell181 and Beam188 elements, and Ref.[41] agree well with each other under different τ_b, which verifies the proposed model.

Based on the different models, the X-, Y-, Z-, and θ-vibration responses of the tip node (see Fig. 3) under the base excitation are shown in Fig. 5, where

It should be noted that the vibration responses of 17 excitation periods are calculated, and each period is divided into 60 time steps. In addition, the Rayleigh damping C is adopted, and the corresponding expressions are listed as follows^[42]:

(16)

where , and . is a range of concerned frequencies, which is closer to the first-two order natural frequencies of the RPICBM, ξ₁ and ξ₂ are the corresponding modal damping ratios of the RPICBM. In this paper, f₁ = 195.56 Hz, f₂= 733.05Hz, and ξ₁=ξ₂=0.03.

Generally speaking, the X-, Y-, and Z-vibration responses agree well with each other (see Figs. 5(a)–5(c)), while the θ-vibration response obtained from the proposed model is closer to that obtained from the Shell181 elements than the Beam188 elements (see Fig. 5(d)), which indicates that there may exist some errors during the process of simulation for the torsional behaviors of the RPICBM when the Beam188 elements are used.

3 Dynamic characteristics of the RPICBM

In this section, the pre-stressed modal^[43], where the rotation speed-induced Coriolis, centrifugal stiffening, and spin softening effects are considered, and the transient dynamic analyses are performed to study the effects of β₀ and γ(L) on the natural frequencies and the vibration responses. The corresponding simulation parameters are listed in Table 1.

3.1 Pre-stressed modal analysis

With the consideration of the rotation speed-induced Coriolis, centrifugal stiffening, and spin softening effects, the effects of β₀ and γ(L) on Y_I, Z_I, and θ_I are analyzed, respectively (see Figs. 6 and 7).

3.1.1 Case 1: the effects of β₀ on Y_I, Z_I, and θ_I

In Fig. 6, it can be seen that when n≠0r·min^-1, Y_I and θ_I always increase while Z_I decreases when β₀ increases. The larger β₀ is, the larger the increment of Y_I and θ_I is (see Figs. 6(a) and 6(c)) and the larger the decrement of Z_I is (see Fig. 6(b)). This is because the spin softening effect decreases in the Y- and θ-directions (see in Eq.(6) and in Eq.(8)) but increases in the Z-direction when β₀ increases (see in Eq.(7)), which leads to the increase in Y_I and θ_I and the decrease in Z_I (see Fig. 6)

3.1.2 Case 2: the effects of γ(L) on Y_I, Z_I and θ_I

Some typical features in Fig. 7 are concluded as follows:

(ⅰ) When γ(L) increases, Y_I and θ_I increase, while Z_I decreases (see Fig. 7). This is because the element structural stiffness increases in the Y- and θ-directions but decreases in the Z-direction when γ(L) increases.

(ⅱ) The larger γ(L) is, the larger the increment of Y_I and θ_I is (see Figs. 7(a) and 7(c)) and the smaller the increment of Z_I is (see Fig. 7(b)). This is because the increasing γ(L) reduces the spin softening effects in the Y- and θ-directions (see in Eq.(6) and in Eq.(8)) while increases the spin softening effect in the Z-direction (see in Eq.(7)). This leads to the increase in Y_I and θ_I and the decrease in Z_I (see Fig. 7).

3.2 Vibration responses in the X-, Y-, Z-, and θ-directions subject to the X-harmonic base excitation

Considering the effect of the X-harmonic base excitation (see Fig. 3), we perform a transient dynamic analysis to obtain the amplitude frequency responses in the X-, Y-, Z-, and θ-directions as well as the constant components in the X- and θ-directions. The corresponding simulation parameter settings can refer to Cases 3 and 4 in Table 1.

3.2.1 Case 3: the effects of β₀ on the vibration responses under the base excitation

The amplitude frequency responses in the X-, Y-, Z-, and θ-directions are shown in Fig. 8. Some main conclusions are made as follows:

(ⅰ) When β₀ = 0° and γ(L)= 0°, the first-order resonance is only embodied in the X- and Y-directions, respectively, (see Figs. 8(a)–8(d)), while the amplitude frequency responses in Figs. 8(c) and 8(d) are zero. This is because there exists the coupling between the X- and Y-directions via the Coriolis effect (see in Eq.(5) and in Eq.(6)), while there is no coupling between the Z- and X- or Y-directions (see in Eq.(5), in Eq.(6), and in Eq.(7)). Moreover, the base excitation load in the Z-direction is zero (see in Eq.(7) and in Eq.(8)). This directly leads to the disappearance of the Z- and θ-vibration responses (see Figs. 8(c) and 8(d)).

(ⅱ) When β₀ = ± 20° and γ(L)=0°, there are couplings among the X-, Y-, Z-, and θ-directions, which leads to the appearance of the resonance peak. The corresponding coupling terms are mentioned in Eq.(1).

(ⅲ) When β₀ = ± 45° and γ(L)=0°, only vibration responses in the X-, Y-, Z-, and θ-directions occur, and there is no resonance peak. This is mainly due to that the excitation frequency n/60 is lower than the first-order flapwise dynamic frequency of the model under β₀ = 45° or -45° and γ(L)=0° (see Fig. 9).

Furthermore, it is worth noting that the constant components in the X- and θ-directions under different β₀ and n are also extracted (see Fig. 10). Due to the non-existence of constant loads in the Y- and Z-directions (see Eqs.(6), (7), (9), and (10)), the constant components in the Y- and Z-directions are zero and omitted here.

In Fig. 10(a), the X-constant components are not affected by β₀, and increase with the increase in n. This is because the element centrifugal force obtained from the proposed model (see in Eq.(5)) is independent of β₀ and increases with the increase in n .

In Fig. 10(b), the θ-constant component is the largest when β₀ = ± 45° and the smallest when β₀ = 0°. Moreover, a positive/negative β₀ leads to a positive/negative θ-constant component (see Figs. 10(b) and 11). This can be easily illustrated via the constant element torque (see Eq.(8), and γ(L) = 0°. When β₀ = ± 45° and 0°, respectively, the magnitude of the constant element torque is the largest and 0 N·m. This leads to the largest and smallest θ-constant components (see Fig. 10(b)).

3.2.2 Case 4: the effects of γ(L) on the vibration responses under the base excitation

Different from the effects of β₀ on the vibration responses under the base excitation (see Subsection 3.2.1), γ(L) (γ(L) ≠ 0°) directly causes couplings among the X-, Y-, Z-, and θ-directions, and further leads to the occurrence of the first-order resonance in all the directions (see Figs. 12 and 13). The corresponding coupling terms can refer to Subsection 3.2.1.

In Fig. 14(a), the X-constant component obtained from the proposed model is not affected by γ(L). This is because the element centrifugal force (see in Eq.(5)) is also independent of γ(L).

In Fig. 14(b), when the summation of two different γ(L) is 0°, the θ-constant component is symmetric about θ = 0 rad. Furthermore, a positive/negative γ(L) can result in a negative/positve θ-constant component (see Figs. 14(b) and 15). The larger γ(L) is, the larger the θ-constant component is. The reasons for these phenomena can be explained as follows:

(Ⅰ) The element torque in the θ-direction is

(see Eq.(8), and β₀ = 0°), and the latter term with a double underline is larger than the former in this study when .

(Ⅱ) A negative/positive γ(L) can determine a positive/negative element torque, thus leading to a positive or negative θ-constant component (see Figs. 5(d), 14(b), and 15).

4 Conclusions

The coupled flapwise-chordwise-axial-torsional dynamic model of the RPICBM is established with the Hamilton principle and the FE method. The proposed model is verified via comparisons with those in the literature. In the FE models, Beam188 and Shell181 elements are used in the ANSYS software. Based on the proposed model, the effects of the setting angle and the pre-twisted angle on the natural frequencies and the base excitation-induced vibration responses of the RPICBM are discussed. Some conclusions are summarized as follows:

(ⅰ) When the setting angle or pre-twisted angle increases, the first-order flapwise and torsional frequencies increase while the first-order chordwise frequency decreases. Moreover, the first-order flapwise, chordwise, and torsional frequencies are affected by the pre-twisted angle rather than the setting angle when the rotation speed is 0 r·min^-1.

(ⅱ) Under the effect of the axial base excitation, there exists the flapwise-axial coupling or even the flapwise-chordwise-axial-torsional coupling of the RPICBM. In addition, a positive or negative setting angle can result in a positive or negative constant component, while a positive or negative pre-twisted angle can result in a negative or positive constant component. When the pre-twisted angle is 0°, and the setting angle is ± 45°, the constant component in the torsional direction reaches the maximum. When both the setting angle and the summation of two pre-twisted angles are 0°, the constant components in the torsional direction under two pre-twisted angles are the same in the magnitude but opposite in the axial directions.

(ⅲ) Within the scope of linear elastic theory, the FE model using Shell181 elements in the ANSYS software and the proposed model may be more suitable for simulating the dynamic characteristics of the RPICBM than the FE model using Beam188 elements in the ANSYS software, especially for the simulation of torsional behaviors.

Appendix A

The derivation of the shape functions in the y^e-bending, z^e-bending, x^e-axial, and x^e-torsional directions are as follows.

(ⅰ) The shape function in the x^e-direction is

(A1)

(ⅱ) The torsional shape function in the θ-direction is

(A2)

(ⅲ) The bending shape function in the x^eo^ey^e-plane is

(A3)

(ⅳ) The bending shape function in the x^eo^ez^e-plane is

(A4)

In Eqs.(A1)-(A4), the expressions of the shape functions N_mi(ξ) (m=u, θ; i=1, 2), N_mk(ξ)(m=v, φ, w, ϕ; k=1, 2, 3, 4) can refer to Ref.[38]. However, there exist minor modifications on the ϕ-displacement shape function, and the modified ϕ -displacement shape function is shown as follows:

(A5)

(A6)

(A7)

(A8)

where

Appendix B

The matrices and vectors related to the rotating Timoshenko beam element are as follows.

The mass matrix M^e is

where

The Coriolis matrix G^e is

where

The structural stiffness matrix K_e^c is

where

The centrifugal stiffening matrix K_c^e is

where

The spin softening matrix K_s^e is

where

The external force vector F^e is

where

Appendix C

For the torsional moment of the inertia of a rectangular beam, the torsional stress function τ(y, z) of an arbitrary rectangular cross-sections is^[39]

(C1)

where h is the length of the short edge of the beam cross-section (see Fig.C1), and τ₁(y, z) is the additional stress correction term considering the effect of h/b.

Furthermore, τ(y, z) should satisfy the Poisson's equation as follows:

(C2)

Substituting Eq.(C1) into Eq.(C2), the following expression can be determined:

(C3)

where τ(y, z) should satisfy the following equations at the boundaries of the rectangular cross-section:

(C4)

Substituting Eq.(C1) into Eq.(C4), the following expressions can be determined:

(C5)

Letting τ₁(y, z)=f(y) f(z) and substituting it into Eq.(C3), the following formula can be determined:

(C6)

where λ represents an arbitrary constant.

According to Eq.(C6), the following two expressions can be determined:

(C7)

The general solutions of Eq.(C7) can be obtained as follows:

(C8)

Referring to the theory of membrane analogy, τ1(y, z) should be an even function with respect to the coordinates y and z. Therefore, τ1(y, z) should be expressed as follows:

(C9)

Substituting Eq.(C9) into the second expression in Eq.(C5), the following equation can be determined:

(C10)

The general solutions of Eq.(C10) can be expressed as follows:

(C11)

Substituting Eq.(C11) into Eq.(C9), τ1(y, z) can be expanded in the following series:

(C12)

Substituting Eq.(C12) into the first expression in Eq.(C5), we can obtain

(C13)

Multiplying both sides of Eq.(C13) by cos(λ_ly) and then integrating it with respect to y over [-b/2, b/2], we have

(C14)

Substituting Eq.(C14) into Eq.(C12), we have

(C15)

Combining Eq.(C1) with Eq.(C15), we can obtain the torsional moment of inertia J_s as follows:

(C16)

where

is a correction factor for J_s. In this paper,

when l=0.