© T. R. B. Trevilato et al., Published by EDP Sciences, 2025

1 Introduction

Currently, there is a burgeoning global interest in optimizing industrial activities to make them more sustainable and environmentally friendly in response to significant climate changes, coupled with global demands for economically viable engineering projects while maintaining safety and comfort. This is further reinforced by the decreasing costs of renewable energy generation, that become economically viable as compared to fossil fuelbased energy generation [18]. This context has encouraged the improvement of the performance of rotating machines, and therefore current engineering rotors’ design should not only prioritize efficiency and reliability but also strive for self-sufficiency and sustainability whenever feasible.

The current industrial landscape also urges engineering disciplines to transcend traditional uses of rotors with well-established linear regimes and explore nonlinear ones. Fixed-base rotors experiencing linear operating regimes have been extensively investigated, and their dynamic behavior is well-documented [13]. However, a broad range of applications involve nonlinear rotors subjected to base excitations, such as aircraft engines. Nonetheless, the number of studies addressing this specific area is quite limited [4,5]. Research endeavors related to onboard rotors commonly focus on seismic excitations and excitations transmitted by the base from neighboring machinery and equipment, as claimed by Dakel et al. [8] and Al-Solihat et al. [1].

In more recent studies, Liu et al. [15] investigated the influence of base excitations on a rotor supported by hydrodynamic bearings. The FE model of the analyzed rotor is composed by a flexible shaft, a rigid disk, and two hydrodynamic bearings. The effects of the amplitude and excitation frequency of the base on the dynamic behavior of the rotor were investigated, and the authors demonstrated that the base velocity should not be disregarded when determining the bearing’s hydrodynamic forces. Furthermore, Del Claro et al. [9] evaluated two models for composite material onboard rotors: the first based on the Classical Beam Theory and the other is based on the Unified Theory. Three experimental test benches were employed to validate these methodologies, resulting a complete evaluation of the dynamic behavior of the rotor.

Regarding nonlinear effects affecting the dynamic behavior of rotors, Dakel et al. [7] delved into the nonlinear dynamic behavior of an onboard rotor supported by hydrodynamic bearings. The authors used a rotor FE model, which is able to incorporate various shafts and disks geometries (symmetrical or asymmetrical components). The considered model encompasses four degrees of freedom six types of base motion (three translations and three rotations) and accounts for the nonlinear characteristics of the oil film of the bearing. The Newmark integrator was employed to solve the time-dependent rotor equations of motion.

Lara-Molina et al. [14] conducted a sensitivity analysis of a flexible rotor affected by uncertain parameters, focusing on their effects on the rotor orbits and Campbell diagram. Barbosa et al. [2] carried out a sensitivity analysis of a flexible rotor made of composite material, as modeled by using the Simplified Homogenized Beam Theory [17].

In this context, the present contribution aims to investigate the nonlinear dynamic behavior of an onboard rotor affected by uncertainties in the bearing’s stiffness. A FE model of the rotor was formulated by considering the dynamic forces resulting from the base excitation applied to the onboard rotor [12]. In this case, nonlinearities are introduced in the rotor by considering that the transverse displacements of the shaft are no longer considered as being small (although still within the elastic material regime). The introduction of these nonlinearities and the solution of the resulting equations of motion follow the approach proposed by Zienkiewicz and Taylor [20]. The proposed model incorporates displacements in the transverse directions as well as rotations around these directions, which have yielded satisfactory results, as will be demonstrated later. Finally, uncertainties are included in the stiffness of the bearings, where the random field is discretized using the Karhunen-Loeve expansion method, thus generating uncertain envelopes for the results in the time domain.

It is worth mentioning that the novelty of the present contribution is the numerical and experimental investigations regarding the dynamic behavior of an onboard rotor in the existence of nonlinearities and uncertainties affecting the dynamic behavior of the system.

Fig. 1 Schematic representation of an onboard rotor.

2 Mathematical model

This section presents the mathematical model of onboard rotors based on work of [10]. Part of the mathematical development is reproduced here, as it is relevant to the progression of this study and contributes to making the article more self-contained.

Taking the schematic representation of an onboard rotor shown in Figure 1 as a reference, wherein three reference axes are employed [10], R₀(x₀, y₀, z₀) denotes the inertial coordinate system, R_S(x_S, y_S, z_S) represents the mobile frame fixed to the rotor base and R(x, y, z) characterizes the mobile frame fixed to the disk.

The rotor orientations are described by using the Euler angles (Fig. 2) [13]. In this case, ψ, θ, and ϕ represent the motion of the reference R relative to the coordinate system R_S, while α, β, and λ stand for the motion of R_S relative to R₀. Figure 2a represents the following reference rotations:

• A rotation ψ around z_S results in the intermediate reference frame R₁ (x₁, y₁, z₁);

• A rotation θ around leads to another intermediate reference frame R₂(x₂, y₂, z₂);

• A rotation ϕ is performed around y with y//y₂.

Similarly, based on Figure 2b, the following can be deduced:

• A rotation α around z₀ results in the intermediate reference frame R₃ (x₃, y₃, z₃);

• A rotation β around x₃ leads to another intermediate reference frame R₄ (x₄, y₄, z₄);

• A rotation λ is performed around y with y//y₃.

Therefore, the angular velocity of R as described in the reference frame R_S is expressed by equation (1). The angular velocity of R_S as described in the reference frame R₀ is represented by equation (2). The determination of the angular velocity of R (disk) as referenced in frame R₀ (ground) is described as given by equation (3).

$${\overrightarrow{\text{Ω}}}_R^{R_S}\text{=}\left\{\begin{array}{c}\dot{\theta }\hspace{0.17em}\text{cos}\hspace{0.17em}\varphi -\dot{\psi }\hspace{0.17em}\text{cos}\hspace{0.17em}\theta \hspace{0.17em}\text{sin}\hspace{0.17em}\varphi \\ \dot{\varphi }+\dot{\psi }\hspace{0.17em}\text{sin}\hspace{0.17em}\theta \\ \dot{\theta }\hspace{0.17em}\text{sin}\hspace{0.17em}\varphi +\dot{\psi }\hspace{0.17em}\text{cos}\hspace{0.17em}\theta \hspace{0.17em}\text{cos}\hspace{0.17em}\varphi \end{array}\right\}$$(1)

$${\overrightarrow{\text{Ω}}}_{R_S}^{R_{\text{0}}}\text{=}\left\{\begin{array}{c}\dot{\beta }\hspace{0.17em}\text{cos}\hspace{0.17em}\gamma -\dot{\alpha }\hspace{0.17em}\text{cos}\hspace{0.17em}\beta \hspace{0.17em}\sin \gamma \\ \dot{\gamma }+\dot{\alpha }\sin \beta \\ \dot{\beta }\hspace{0.17em}\text{sin}\hspace{0.17em}\gamma +\dot{\alpha }\hspace{0.17em}\cos \hspace{0.17em}\beta \hspace{0.17em}\cos \gamma \end{array}\right\}\text{=}\left\{\begin{array}{c}{\dot{\alpha }}_S\\ {\dot{\beta }}_S\\ {\dot{\gamma }}_S\end{array}\right\}$$(2)

$${\overrightarrow{\text{Ω}}}_R^{R_S}={\overrightarrow{\text{Ω}}}_{R_S}^{R_{\text{0}}}+{\overrightarrow{\text{Ω}}}_R^{R_S}\text{=}\left\{\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right\}$$(3)

where $\left({}^{\cdot }\right)$ denotes derivative with respect to time. The spatial position of the rotor base (vector $\overline{OA}$; Fig. 1) concerning the reference frame R₀ is described by the coordinates x_a, y_a, and z_a, as given by equation (4). Conversely, the coordinates X, Y, and Z represent the position of this identical vector concerning the reference frame R_S (Eq. (5)).

$${\overline{OA}}^{R0}=\left\{\begin{array}{c}{x}_A\\ y_A\\ z_A\end{array}\right\}$$(4)

$$\begin{array}{l}{\overline{OA}}^{R_S}\text{=}\left\{\begin{array}{c}X\\ Y\\ Z\end{array}\right\}\\ =\left\{\begin{array}{c}\left(x_A\cos \alpha +y_A\sin \alpha \right)\cos \gamma -\left[z_A\cos \beta +\left(x_A\sin \alpha -y_A\cos \alpha \right)\sin \beta \right]\sin \gamma \\ z_A\sin \beta -\left(x_A\sin \alpha -y_A\cos \alpha \right)\cos \beta \\ \left(x_A\cos \alpha +y_A\sin \alpha \right)\sin \gamma -\left[z_A\cos \beta +\left(x_A\sin \alpha -y_A\cos \alpha \right)\sin \beta \right]\cos \gamma \end{array}\right\}\end{array}$$(5)

The position of the point C (vector $\overline{AC}$; Fig. 1) with respect to the reference frame R_S is presented by equation (6). Equation (7) presents the vector $\overline{OC}$ described in the reference frame R_S.

$${\overline{AC}}^{R_{\text{S}}}\text{=}\left\{\begin{array}{c}u\text{(}y\text{,}t\text{)}\\ y\\ w\text{(}y\text{,}t\text{)}\end{array}\right\}$$(6)

$${\overline{OC}}^{R_S}={\overline{OA}}^{R_S}+{\overline{AC}}^{R_S}=\left\{\begin{array}{c}X+u\text{(}y\text{,}t\text{)}\\ Y+y\\ Z+w\text{(}y\text{,}t\text{)}\end{array}\right\}$$(7)

herein u and w represent the transverse displacements of the shaft in R_S along with the x_S and z_S axes, respectively. Additionally, y denotes the spatial coordinate of the disk along y_S.

Fig. 2 Transformation of coordinates (Adapted from Lalanne and Ferraris [13]).

2.1 Kinetic energy

The kinetic energies of the symmetrical disk (T_D), the unbalance mass (T_u), and the symmetrical shaft (T_S) are provided by equations (8), (9), and (10), respectively.

$$T_D\text{=}\frac{1}{\text{2}}{m}_D{\left({\overrightarrow{V_D}}^{R_0}\right)}^{\text{2}}+\frac{1}{\text{2}}{\overrightarrow{\text{Ω}}}_R^{R_0}{I}_D{\overrightarrow{\text{Ω}}}_R^{R_0}$$(8)

$$T_u\text{=}\frac{1}{\text{2}}{m}_u{\left({\overrightarrow{V_D}}^{R_0}\right)}^2$$(9)

$$\begin{array}{l}{T}_S\text{=}\frac{1}{\text{2}}\rho S{\displaystyle {\int }_0^L\text{(}{\dot{u}}_E^2+{\dot{v}}_E^2+{\dot{w}}_E^2\text{)d}y}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{+}\frac{1}{\text{2}}\rho {\displaystyle {\int }_0^L\text{(}{I}_x{\omega }_x^2+I_y{\omega }_y^2+I_z{\omega }_z^2\text{)d}y}\end{array}$$(10)

where I_x, I_y, and I_z are the area moments of inertia of the shaft associated with the axes x, y, and z, correspondingly. Meanwhile, m_D is the mass of the disk, ρ denotes the material density, S is the cross-sectional area of the shaft, L is the shaft length, and m_u represents the unbalance mass.

The translation velocity V_C of the disk and the velocity V_D of the unbalance mass, both relative to reference frame R_S, are given by equations (11) and (12), respectively.

$$\begin{array}{l}{V}_C^{R_S}\text{=}\frac{\text{d}}{\text{d}t}\left({\overline{OC}}^{R_S}\right)+{\overrightarrow{\text{Ω}}}_{R_S}^{R_0}\times {\overline{OC}}^{R_{\text{S}}}\\ ={\left\{\begin{array}{c}\text{(}\dot{X}+\dot{u}+{\dot{\beta }}_S\text{(}Z+w\text{)}-{\dot{\gamma }}_S\left(Y+y\right)\\ \dot{Y}+{\dot{\gamma }}_S\text{(}X+u\text{)}-{\dot{\alpha }}_S\text{(}Z+w\text{)}\\ \dot{Z}+\dot{w}+{\dot{\alpha }}_S\left(Y+y\right)-{\beta }_S\left(X+u\right)\end{array}\right\}}_{R_s}\text{=}{\left\{\begin{array}{c}{\dot{u}}_C\\ {\dot{v}}_C\\ {\dot{w}}_C\end{array}\right\}}_{R_s}\end{array}$$(11)

$$\begin{array}{l}{V}_D^{R_S}\text{=}\frac{\text{d}}{\text{d}t}\left({\overline{OC}}^{R_S}\right)+{\overrightarrow{\text{Ω}}}_{R_S}^{R_0}\times {\overline{OC}}^{R_{\text{S}}}\\ \hspace{0.17em}\hspace{0.17em}\text{=}{\left\{\begin{array}{c}\dot{X}+\dot{u}+\text{dΩ}\hspace{0.17em}\cos \text{Ω}t+{\dot{\beta }}_S\text{(}Z+w+d\text{Ω}\cos \hspace{0.17em}\text{Ω}t\text{)}-{\dot{\gamma }}_S\left(Y+y\right)\\ \dot{Y}+{\dot{\gamma }}_S\text{(}X+u+\text{dΩ}\hspace{0.17em}\sin \text{Ω}t\text{)}-{\dot{\alpha }}_S\text{(}Z+w+\text{dΩ}\cos \text{Ω}t\text{)}\\ \dot{Z}+\dot{w}-\text{dΩ}\sin \text{Ω}t+{\dot{\alpha }}_S\left(Y+y\right)-{\dot{\beta }}_S\left(X+u-\text{dΩ}\sin \text{Ω}t\right)\end{array}\right\}}_{R_s}\end{array}$$(12)

where t stands for the time and d represents the distance of the unbalance mass m_u from the geometric center of the shaft at the position C and the point D represents the location of the unbalance mass, with the vector $\overline{OD}$ relative to the reference frame R_S given by:

$${\overline{OD}}^{R_S}={\overline{OC}}^{R_S}+{\overline{CD}}^{R_S}=\left\{\begin{array}{c}X+u\text{(}y,t\text{)}+\text{d}\cos \hspace{0.17em}\text{Ω}t\\ Y+y\\ Z+w\text{(}y,t\text{)}+\text{d}\sin \text{Ω}t\end{array}\right\}$$(13)

The mass moment of inertia tensor of the disk, denoted as I_D (see Eq. (14)), comprises I_Dx, I_Dy, and I_Dz components, representing the mass moments of inertia associated with the x, y, and z directions, respectively.

$$I_D\text{=}\left\{\begin{array}{ccc}{I}_{Dx}& 0& 0\\ 0& I_{Dy}& 0\\ 0& 0& I_{Dz}\end{array}\right\}$$(14)

2.2 Strain energy

The strain energy U of the shaft remains constant along the motion of the rotor’s base [19]. This is associated to the fact that the deformation of the shaft results from the constraints imposed by the bearings. Therefore, for the formulation of the strain energy it is imperative to describe the shaft displacement field, as given by equation (15).

$$\widehat{u}=\left\{\begin{array}{c}{u}_{\text{e}}\\ v_{\text{e}}\\ w_{\text{e}}\end{array}\right\}\text{=}\left\{\begin{array}{c}{u}_0\text{(}y,t\text{)}\\ 0\\ w_0\text{(}y,t\text{)}\end{array}\right\}+\left\{\begin{array}{c}0\\ -x\frac{\partial u_e}{\partial y}-z\frac{\partial w_e}{\partial y}\\ 0\end{array}\right\}$$(15)

where u_e, v_e, and w_e stand for the displacements in the x, y, and z directions, respectively, and u₀ and w₀ represent the displacements of the shaft mid-line along with the x and z directions, respectively.

Applying the displacement field to the Green-Lagrange strain tensor,

$$e\left(\widehat{u}\right)=\frac{1}{\text{2}}\text{(}\nabla \widehat{u}+\nabla {\widehat{u}}^T\text{)}+\frac{1}{\text{2}}\nabla \widehat{u}\nabla {\widehat{u}}^T$$(16)

in which ∇ is the gradient operator.

In scenarios characterized by small displacements the term 1/2 (∇u∇u^T) can be neglected. However, in cases involving significant displacements, this term cannot be disregarded [12,20]. Thus, the strain field retaining the significant terms is given by:

$$\begin{array}{l}\text{e}=\left\{\begin{array}{c}{e}_{yy}\\ e_{xy}\\ e_{yz}\end{array}\right\}=z\left\{\begin{array}{c}-w_{\text{0,}yy}\\ 0\\ 0\end{array}\right\}+x\left\{\begin{array}{c}-u_{\text{0,}yy}\\ 0\\ 0\end{array}\right\}\\ +z\left\{\begin{array}{c}0\\ u_{\text{0,}y}{w}_{\text{0,}yy}\\ w_{\text{0,}y}{w}_{\text{0,}yy}\end{array}\right\}+x\left\{\begin{array}{c}0\\ u_{\text{0,}y}{u}_{\text{0,}yy}\\ w_{\text{0,}y}{u}_{\text{0,}yy}\end{array}\right\}+\dots \\ \hspace{1em}\text{…}+\frac{1}{\text{2}}{z}^2\left\{\begin{array}{c}{w}_{\text{0,}yy}^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}\\ 0\\ 0\end{array}\right\}+\frac{1}{\text{2}}{x}^2\left\{\begin{array}{c}{u}_{0,yy}^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}\\ 0\\ 0\end{array}\right\}\\ \hspace{1em}\text{+}\frac{1}{\text{2}}\left\{\begin{array}{c}{w}_{\text{0,}y}^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}+w_{\text{0,}y}^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}\\ 0\\ 0\end{array}\right\}\end{array}$$(17)

where {(),_y = ∂ ()/∂y}.

Therefore, the strain energy of the system can be written as follows:

$$U=\frac{1}{\text{2}}{{\displaystyle \int }}^{\text{}}\left(E{e}_{yy}^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}+G\left(e_{xy}^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}+e_{yz}^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}\right)\right)dV$$(18)

where E and G are the Young and shear modulus, respectively.

2.3 Shaft FE model

The FE model of the shaft assumes a circular cross-section with constant diameter along the shaft length L, as presented in Figure 3. In this case, a beam theory is used to formulate a FE element with two nodes and four degrees of freedom each node: two displacements and two rotations (displacements u₁, u₂, w₁, and w₂; and rotations θ₁, θ₂, ψ₁, and ψ₂). The vector of nodal displacements for the shaft is given by equation (19).

$$q={\left[u_1\hspace{1em}{w}_1\hspace{1em}{\theta }_1\hspace{1em}{\psi }_1\hspace{1em}{u}_2\hspace{1em}{w}_2\hspace{1em}{\theta }_2\hspace{1em}{\psi }_2\right]}^T$$(19)

Consequently, the equations of motion representing the dynamic behavior of the rotor is written as:

$$\begin{array}{l}\left[M_D+M_E\left]\ddot{q}+\hspace{0.17em}\right[D_D+D_D^{\text{*}}+D_E+D_b\right]\dot{q}+\dots \\ \dots +\left[K_D^{\text{*}}+K_E^{\text{*}}+K_E+K_{nl}+K_b\right]q=W+F_u+F_D^{\text{*}}+F_E^{\text{*}}\end{array}$$(20)

in which M_D and D_D represent the classical mass and gyroscopic matrices of the disk, respectively. The matrices M_E and D_E are the classic matrices of mass and gyroscopic of the shaft, respectively, while K_E and K_nl are the linear and nonlinear (representing the bending coupling between the x_S and z_S directions) stiffness matrices of the shaft, respectively. The matrices $D_D^{*},K_D^{*}$, and $K_E^{*}$ are associated with the motion of the rotor base, as well as the force vectors $F_D^{*}$ and $F_E^{*}$. The matrices K_b and D_b stand for the stiffness and damping due to the bearings, W represents the weight of the rotating components, and F_u is the unbalance forces. A comprehensive description of the mentioned FE matrices can be found in the “Appendix” of the present contribution.

Fig. 3 Degrees of freedom associated with the finite element of the shaft [2].

2.4 Stochastic FE model

In the context of this study, the stiffness coefficients of the bearings were treated as uncertain parameters following a Gaussian random distribution [3,11]. To compute the associated stochastic stiffness matrices, the proposed model employs the Karhunen-Loève (KL) expansion method and the Latin Hypercube sampling to obtain a stochastic distribution of the mechanical properties. Given a stochastic process H(x, z, Θ),

$$\overline{H}\left(x,z,\text{Θ}\right)=E\left[H\left(x,z,\text{Θ}\right)\right]$$(21)

$$\begin{array}{l}C\left[H\text{(}{x}_1,z_1,\text{Θ)},H\text{(}{x}_2,z_2,\text{Θ)}\right]=\dots \\ \text{…}=E\left[H\text{(}{x}_1,z_1,\text{Θ)}-\overline{H}\text{(}{x}_1,z_1\text{)}\left]E\right[H\text{(}{x}_2,z_2,\text{Θ)}-\overline{H}\text{(}{x}_2,z_2\text{)}\right]\end{array}$$(22)

where Θ is the variable associated to a realization (or sample) of the stochastic process, E [•] is the expectation operator, and C [•] stands for the covariance. Consequently, assuming that H(x, z, Θ) is a homogeneous Gaussian process with a symmetric and positive-definite covariance within the physical domain, the finite statistical approximation of this process of order m is expressed as follows [11]:

$$H\text{(}x,z,\text{Θ)}\approx \overline{H}\text{(}x,z\text{)}+\underset{r=1}{\overset{m}{{{\displaystyle \sum }}^{\text{}}}}\sqrt{{\lambda }_r}{f}_r\text{(}x,z\text{)}{\xi }_r\text{(Θ)}$$(23)

where ξ_r (Θ) represents uncorrelated Gaussian variables. The terms λ_r = λ_iλ_j and f_r (x, z) = f_i (x) f_j (y) correspond to the eigenvalues and eigenfunctions defined through the covariance function, respectively.

The eigenvalues and eigenfunctions are determined in terms of the roots w_r of the following odd and even transcendental equations.

For odd i and j, with i ≥ 1 and j ≥ 1:

$$\begin{array}{l}{\lambda }_i=\frac{2c_1}{w_i{}^2+c_1{}^2},\hspace{1em}{f}_i\text{(}x\text{)}={\alpha }_i\cos \text{(}{w}_{i}x\text{)}\\ {\lambda }_i=\frac{2c_2}{w_j{}^2+c_2{}^2},\hspace{1em}{f}_j\text{(}z\text{)}={\alpha }_j\cos \text{(}{w}_{j}z\text{)}\end{array}$$(24)

with:

$$\begin{array}{l}{\alpha }_i=\frac{1}{\sqrt{a+\sin \text{(}{w}_{i}a\text{)/2}{w}_i}},\hspace{1em}{\alpha }_j\text{=}\frac{1}{\sqrt{b+\sin \text{(}{w}_{j}b\text{)/2}{w}_j}}\\ c_1=\frac{1}{L_{cor,x}},\hspace{1em}{c}_2=\frac{1}{L_{cor,z}}\end{array}$$(25)

Here, L_cor,x and L_cor,z represent the correlation lengths in the x and z directions within domains [−a, a] and [−a, a] where a and b denote the dimensions of the finite element used in the x_S and z_S directions, respectively. In the specific case of this work, both a and b are equal to the radius of the shaft.

The transcendental equations and the domain within which they are defined are obtained as follows:

$$\begin{array}{l}\text{Equation}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{:}{c}_1-w_i\tan \text{(}{w}_{i}a\text{)}=\text{0;}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{c}_2-w_j\tan \text{(}{w}_{j}b\text{)}=0\\ \text{Domain}\hspace{0.17em}\hspace{0.17em}\text{:}\left[\text{(}i-\text{1)}\pi \text{/}a,\text{(}i-\text{1/2)}\pi \text{/}a\right]\hspace{1em}\hspace{1em}\left[\text{(}j-\text{1)}\pi \text{/}a,\text{(}j-\text{1/2)}\pi \text{/}a\right]\end{array}$$(26)

For i and j even, with i ≥ 2 and j ≥ 2:

$$\begin{array}{l}{\lambda }_i\text{=}\frac{2c_1}{w_i^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}+c_1^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}},\hspace{1em}{f}_i\text{(}x\text{)}={\alpha }_i\sin \text{(}{w}_{i}x\text{)}\\ {\lambda }_i\text{=}\frac{2c_2}{w_j^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}+c_2^{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}},\hspace{1em}{f}_j\text{(}x\text{)}={\alpha }_j\sin \text{(}{w}_{j}x\text{)}\end{array}$$(27)

with:

$${\alpha }_i=\frac{1}{\sqrt{a-\sin \left(w_{i}a\right)/2w_i}},\hspace{1em}{\alpha }_j=\frac{1}{\sqrt{b-\sin \left(w_{j}b\right)/2w_j}}$$(28)

The transcendental equations and the domain in which they are defined are given by:

$$\begin{array}{l}\text{Equation}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{:}{c}_1+w_i\tan \text{(}{w}_{i}a\text{)}=\text{0,}\hspace{1em}\hspace{1em}{c}_2-w_j\tan \text{(}{w}_{j}b\text{)}=0\\ \text{Domain}\hspace{0.17em}\hspace{0.17em}\text{:}\left[\text{(}i-\text{1/2)}\pi \text{/}a,\hspace{0.17em}\hspace{0.17em}i\pi \text{/}a\right]\hspace{1em}\hspace{1em}\left[\text{(}j-\text{1/2)}\pi \text{/}a,j\pi \text{/}a\right]\end{array}$$(29)

The KL expansion has been employed for modeling the elementary random stiffness matrix of the bearings as given by:

$$\text{Θ}{K}_b\text{(Θ)}={\overline{K}}_b\text{+}\underset{r=1}{\overset{m}{{{\displaystyle \sum }}^{\text{}}}}{{\text{[}}_{st}{K}_b\text{]}}^{\text{(}k\text{)}}{\xi }_r^k\text{(Θ)}$$(30)

where ${\overline{K}}_b$ is the mean matrix calculated according to equation (20), while the stochastic matrix is defined as follows:

$${}_{st}{K}_b={\alpha }_i{\alpha }_j\sqrt{{\lambda }_i}\sqrt{{\lambda }_j}\underset{x=a}{\overset{x=0}{{{\displaystyle \int }}^{\text{}}}}\underset{z=b}{\overset{z=0}{{{\displaystyle \int }}^{\text{}}}}{f}_i\text{(}x\text{)}{f}_j\text{(}z\text{)}\left\{\begin{array}{cc}{k}_{xx}& k_{xz}\\ k_{zx}& k_{zz}\end{array}\right\}\text{d}x\text{d}z$$(31)

in which k_xx, k_xz, k_zx, and k_zz are the stiffness coefficients of the bearings along the x and z directions.

3 Rotor test rig

Figure 4a shows the test rig used in the study of the onboard rotor. The rotor consists of a 550 mm long and 10 mm diameter flexible steel shaft (FE model composed of 21 beam elements; density ρ = 7800 kg/m³, Young’s modulus E = 190 GPa, and Poisson ratio v = 0.3), a rigid disk D (node 10, 2.3 kg mass) made of steel with 100 mm diameter and 40mm thickness (ρ = 7842 kg/m³), and two self-alignment ball bearings (B₁ and B₂, located at nodes 3 and 20, respectively). Displacement sensors (Meggitt TQ-412) were used to collect the rotor vibration responses at the disk position (orthogonally placed at node 10; sensors S_10xs and S_10zs). The system is driven by a DC electric motor. Measurements were taken using an Agilent® analyzer (model 35670A) with a range of 0 to 250 Hz in steps of 0.25 Hz. The rotational speed was controlled through a power source and verified by an encoder.

A model parameter updating procedure was performed to obtain the unknown parameters to be used in the FE model of the rotor, including the stiffness and damping coefficients of the bearings, the stiffness due to coupling between the electric motor and the shaft (K_rotx and K_rotz, around the orthogonal directions of the node 1), and the modal damping coefficients for the first two vibration modes.

The updating process was performed through the solution of a typical inverse problem based on simulated and experimental Frequency Response Functions (FRFs) of the rotor [6]. For this aim, the Differential Evolution optimization approach was applied [16].The objective function adopted in this case is presented by equation (32), in which only regions in the vicinity of the peaks associated with the natural frequencies were taken into consideration. The minimization of the objective function was carried out with an initial population of 100 individuals, and this process was repeated 10 times to determine the global minimum solution of the problem.

$$Objective={\displaystyle \sum _{i=1}^n\frac{\Vert FR{F}_{\text{exp,}i}-FR{F}_{num\text{,}i}\Vert }{\Vert FR{F}_{\text{exp,}i}\Vert }}$$(32)

where n is the number of FRFs used in the procedure, FRF_exp represents the experimental data, and FRF_num corresponds to the numerical results obtained by the FE model of the rotating system.

The FRFs were acquired with the rotor at rest by applying impact forces along the x_S and z_S directions of the disk. The signal was measured by a proximity sensor positioned in the direction of the applied force, as shown in Figure 5, resulting two FRFs. Table 1 presents the parameters determined in the end of the optimization process, and Figure 5 illustrates the comparisons between the experimental FRFs and those obtained by the updated FE model. Note that the FE model is representative, particularly in the frequency range that encompasses the first two vibration modes of the rotor. The critical speeds of the rotor are approximately 1620 RPM (BW) and 1640 RPM (FW), as shown in the Campbell diagram presented in Figure 6.

Fig. 4 Rotor test rig and FE model.

4 Deterministic analyses

This section is dedicated to evaluating the dynamic behavior of the rotor under different base excitations.

In this case, the excitations were performed by using the Dongling Vibration® electrodynamic shaker (model GT1000) presented in Figure 7. In the conducted tests, the excitations were applied in the xS direction with the rotor operating at 1200 RPM (excitations measured by using an accelerometer).

Initially, the rotor was subjected to impacts with two different amplitudes: 1.5 m/s² and 25 m/s² within a 5 ms interval (Fig. 8). The associated numerical and experimental vibration responses of the rotor are presented in Figures 9 and 10, respectively. The vibration responses of the linear and nonlinear FE models are presented in Figures 9 and 10 for comparison purposes. Note that the 1.5 m/s² excitation amplitude is insufficient to reveal the nonlinear existing in the rotor (see Fig. 9a). Consequently, the responses corresponding to the linear and nonlinear FE models were nearly identical. Furthermore, it can be concluded that the FE model adequately represents the dynamic behavior of the rotor in both directions for this excitation amplitude.

Figure 10 illustrates the vibration responses of the rotor obtained for a base excitation of 25 m/s² amplitude. Note that the nonlinear effects became apparent, and the nonlinear FE model demonstrated to be more representative than the linear model. Figure 11 demonstrates the higher representativeness of the nonlinear model can be supported by the Discrete Fourier Transform (DFT) of the vibration responses, as presented in Figure 10. Note that the linear model is not able to represent satisfactorily the damped natural frequency of the rotor at, approximately, 1620 RPM (27 Hz). The unbalance response at 1200 RPM (20 Hz) is adequately represented by both linear and nonlinear FE models. The better response of the nonlinear model is primarily due to its better representation of stiffness increase, typical of nonlinearity of large displacements [20]. It is worth noting that, despite the displacements being smaller than the shaft diameter, It has already entered the large displacement regime, as evidenced by the responses.

Equation (33) presents the sinusoidal excitation applied along with the x_S direction of the rotor.

$${\ddot{x}}_S=\text{Γ}\sin \left(n\frac{2\pi {\text{Ω}}_{cr}}{60}t\right)$$(33)

where Ω_cr = 1560 RPM is a frequency close to the first BW and FW critical speeds of the rotor, n = 1/2 (or n = 1/3) is a constant used to produce submultiple excitations of Ω_cr, and Γ defines the excitation amplitude (15 m/s² and 30 m/s² for n = 1/2 and n = 1/3, respectively). These values were determined based on the limitations of the electrodynamic shaker. Figure 12 presents the applied excitations.

Figure 13 illustrates the vibration responses of the rotor under the applied sinusoidal base excitations. Figures 13a and 13c demonstrate the ability of the nonlinear FE model to represent the dynamic behavior of the rotor in the direction that coincides with the direction of the applied excitation (x_S direction). However, both linear and nonlinear models were not able to reproduce the experimental vibration responses of the rotor along the z_S direction for n = 1/3 and n = 1/2 (see Figs. 13b and 13d). This difficult is associated with the model adopted to represent the coupling between the shaft and the electric motor. The linear stiffnesses K_rotx and K_rotz considered in the formulation are not representative enough in the context of the applied base excitation (see Fig. 12).

Figure 14 illustrates the numerical and experimental DFTs associated with the time domain vibration responses shown in Figure 13. Note that the nonlinear model was able to represent the experimental x_S vibration responses of the rotor for the excitation frequencies corresponding to n = 1/3 and n = 1/2 and at Ω_cr = 1560 RPM (see Figs. 14a and 14c, respectively). The lack of representativeness for both linear and nonlinear models regarding the experimental vibration responses of the rotor along the z_S direction can be observed in Figures 14b and 14d, in which the higher differences between the numerical and experimental curves are observed in the frequency range 1400 to 2000 RPM. However, the nonlinear model can be considered as being better to describe the vibration responses of the rotor (see the performance of the nonlinear model at Ω_cr = 1560 RPM in Fig. 14d).

The Mean Squared Errors (MSE values) of equation (34) considering the numerical and experimental vibration responses were determined for both linear and nonlinear cases for n = 1/2 and n = 1/3.

$$MSE=\frac{1}{N_x}{\displaystyle \sum _{i=1}^{N_x}{\text{(}{\widehat{y}}_i-y_i\text{)}}^2}$$(34)

where N_x represents the number of data points evaluated in the time interval 4s to 4.5 sec, ŷ_i denotes the experimental vibration amplitude and y_i stands for the corresponding numerical responses determined by using either the linear or the nonlinear FE models. Table 2 presents the obtained MSE values and the corresponding standard deviations.

Table 2 reveals that for both cases (n = 1/3 and n = 1/2), the MSE values associated with the nonlinear model are smaller than the ones determined by using the linear model. For the case of n = 1/3, the errors of the nonlinear model are, approximately, 3.6 and 2 times smaller than those of the linear model for the x_S and z_S directions, respectively. In addition, considering n = 1/2, the errors of the nonlinear model were 2.6 and 1.8 times smaller than the errors of the linear model for the x_S and z_S directions, respectively.

Fig. 5 Numerical and experimental FRFs of the rotor from node 10.

Fig. 6 Campbell diagram of the updated FE model.

Fig. 7 Rotating system mounted on the electrodynamic shaker (control accelerometer used in the closed loop vibration control strategy of the shaker).

Fig. 8 Acceleration applied to the base of the rotor along the xS direction.

Fig. 9 Vibration response considering the 1.5 m/s2 impact excitation.

Fig. 10 Vibration response considering the 25 m/s2 impact excitation.

Fig. 11 DFTs of the impact response of 25 m/s2 amplitude.

Fig. 12 Sinusoidal accelerations applied to the rotor base.

Fig. 13 Vibration responses of the rotor under sinusoidal base excitations.

Fig. 14 Numerical and experimental DFTs associated with the vibration responses as presented in Figure 13.

5 Stochastic analyses

However, it is necessary to ascertain whether the MSE values obtained using the nonlinear FE model (smaller MSE values) are statistically relevant. Hypothesis tests were performed for the four scenarios described in Table 2 (n = 1/2 and n = 1/3 for the x_S and z_S directions). The null hypothesis states that the MSE value of the linear model is smaller than or equal to the error obtained by using the nonlinear model, i.e., MSE_L ≤ MSE_NL. Consequently, the alternative hypothesis is MSE_L > MSE_NL. In all test cases, the null hypothesis is rejected in favor of the alternative hypothesis, accepted at a statistical significance level of, approximately, 15%. These results indicate that the nonlinear model leads to a significantly lower error than the linear model with statistical relevance.

Next, a stochastic model was introduced incorporating a 10% level of uncertainty in the bearing stiffness through the so-called Karhunen-Loève expansion (see Section 2.4) (refers to 10% chance in the standard deviation on the related Gaussian variables). In this case, an impulse with 25 m/s² amplitude and a sinusoidal excitation considering n = 1/2 (see details in Eq. (33)) were applied to the rotor base, separately. A total of 200 samples were generated for the analysis of stochastic responses. For both cases the rotor operates at 1200 RPM. A convergence analysis was performed, by using mean square deviation, normalized by the average of all errors, as depicted in Figure 15, demonstrating that the convergence is achieved well before reaching 200 samples, but for greater robustness of the result, 200 samples were maintained.

Figure 16 presents the obtained uncertainty envelopes. As mentioned, only the vibration responses determined by using the nonlinear FE model are evaluated since it reveals to be more adapted to represent the experimental data (see Tab. 2). Concerning the impulse responses (Figs. 16a and 16b), it can be observed that the experimental vibration responses along the xS direction are included in the corresponding uncertainty envelopes. Additionally, the envelopes are narrower immediately after the impulse application and become wider as the signal is attenuated along the time. The uncertain envelopes are not able to represent the experimental vibration responses of the rotor measured along with the z_S direction.

The envelopes associated with the sinusoidal excitation (see Figs. 16c and 16d) is narrower for the vibration responses obtained along with the force direction (x_S direction). Additionally, note that the uncertain envelopes could represent better the vibration responses of the rotor along the x_S direction as compared with the results obtained along the z_S direction. It is also important to highlight that, despite the bearing stiffness being considerably lower than the shaft stiffness, the responses of the nonlinear system were particularly sensitive to uncertainties in this parameter.

Fig. 15 Convergence analyses performed for both base excitation conditions.

Fig. 16 Uncertain envelopes obtained by using the nonlinear FE model under two different excitation conditions.

6 Conclusions

In the present study, a stochastic nonlinear onboard rotor FE model was formulated and experimentally validated. The dynamic analysis of the considered system demonstrates the considerable error margin associated with a linear model when nonlinearity due to large displacements should have been taken into account. In the direction of the applied excitatioArn (x_S direction), the linear model overestimates the vibration response amplitude of the rotor. However, in the z_S direction the linear model fails to capture the coupling of directions due to large displacements, resulting an underestimation of the actual vibration response amplitude. The linear model overestimates the sensitivity to the uncertainties applied to the bearing stiffness. The nonlinear model proposed in this study effectively represents the rotor responses along the x_S direction and shows an improvement in the representation of the vibration responses along the z_S direction as compared with the linear FE model.

It is worth mentioning that some model refinements are scheduled to be introduced in future work, such as the inclusion of shear effects and axial displacement. Finally,uncertainty analysis with respect to various parameters of the FE model by considering the rotor operating under different conditions is also scheduled.

Acknowledgements

The authors are thankful the Brazilian research agencies CNPq (Proc. Nb. 406148/2022-8), FAPEMIG (APQ-040034-24), and CAPES through the INCT-EIE for the financial support provided for this research effort. The authors also would like to thank the companies Foz do Chapecó, Baesa, Enercan and Ceran for technical and financial support, through the Research and Development project PD-02949-3007/2022 with fundings from ANEEL’s R&D program.

Funding

The Article Processing Charges for this article are taken in charge by the French Association of Mechanics (AFM).

Conflicts of interest

The authors declare that they have no conflict of interest.

Data availability statement

The data that support the findings of this study are available on request from the corresponding author, Dr. Aldemir Ap Cavallini Jr (aacjunior@ufu.br).

Author contribution statement

All authors participated and contributed to all parts, from the design and implementation of the research and experimental part, to the analysis of the results and to the writing of the manuscript. All authors agreed to the published version of the manuscript.

Ethical approval

Not applicable.

Appendix A Numerical integration

Equation (20) is numerically integrated by using a variant of the Newmark method. For this aim, the velocity ${\dot{q}}_n$ and acceleration ${\ddot{q}}_n$ vectors are given by [17]:

$$\begin{array}{l}{\ddot{q}}_n\text{=}\frac{1}{{\lambda }_1\text{Δ}{t}^2}\left(q_n-q_{n-1}\right)-\frac{1}{{\lambda }_1\text{Δ}t}{\dot{q}}_{n-1}+\left(1-\frac{1}{2{\lambda }_1}\right){\ddot{q}}_{n-1}\\ {\dot{q}}_n={\dot{q}}_{n-1}+\left(1+{\lambda }_2\right){\ddot{q}}_{n-1}+{\lambda }_2\text{Δ}t{\ddot{q}}_n\end{array}$$(A.1)

where the subscript n indicates that the vector or matrix is calculated at time t_n, and λ₁ and λ₂ are parameters of the Newmark method. Applying equation (A.1) to the equation of motion of the rotor (Eq. (20)), the following equation is obtained:

$$\widehat{K}\left(q_n\right)q_n={\widehat{F}}_n$$(A.2)

where:

$$\begin{array}{l}\widehat{K}\left(q_n\right)=K^{Total}\left(q_n\right)+a_{0}M+a_1{D}^{Total}\\ {\widehat{F}}_n=F_n^{Total}+M{A}_{n-1}+D^{Total}{B}_{n-1}\\ A_{n-1}=a_{0}q{X}_{n-1}+a_2{\dot{q}}_{n-1}+a_3{\ddot{q}}_{n-1}\\ B_{n-1}=a_{1}q{X}_{n-1}+a_4{\dot{q}}_{n-1}+a_5{\ddot{q}}_{n-1}\end{array}$$

with:

$$\begin{array}{l}{a}_0=\frac{1}{{\lambda }_1{\left(\text{Δ}t\right)}^2},\hspace{1em}{a}_1=\frac{{\lambda }_2}{{\lambda }_1\text{Δ}t},\hspace{1em}{a}_2=\frac{1}{{\lambda }_1\text{Δ}t},\\ a_3=\frac{1}{2{\lambda }_1}-1,\hspace{1em}{a}_4=\frac{{\lambda }_2}{{\lambda }_1}-1\hspace{1em}{a}_5=\text{Δ}t\left(\frac{{\lambda }_2}{2{\lambda }_1}-1\right)\end{array}$$

Equation (A.2) is a system of nonlinear algebraic equations at the instant t_n, which can be solved iteratively by using the Newton-Raphson method. This method is based on the expansion of the nonlinear algebraic equation into Taylor series around the known solution. In the formulation used in this work, it is assumed that to obtain the solution at iteration r, the solution at iteration r − 1 is known. Thus, a residual ϕ (q_n) is defined:

$$\varphi \left(q_n\right)\equiv \widehat{K}\left(q_n\right)q_n-{\widehat{F}}_n=0$$(A.3)

Expanding ϕ (q_n) in Taylor series:

$$\varphi \text{(}{q}_n\text{)}=\varphi \left(q_n^{r-1}\right)+\left[\frac{\partial {\varphi }_n^{r-1}}{\partial q_n^{r-1}}\right]\left(q_n^r-q_n^{r-1}\right)+\text{O}\left(q^2\right)$$(A.4)

where O q² represents higher order terms. Using equation (A.3) and (A.4) and by considering the higher order terms as negligible:

$${}^{t}K\left(q_n^{r-1}\right)\text{Δ}{q}^r+{\varphi }_n^{r-1}=0$$(A.5)

where ^tK is known as tangent matrix, or tangent stiffness matrix. The displacement increment can be given by:

$$q_n^r=q_n^{r-1}+\text{Δ}{q}^r$$(A.6)

The calculation of the ^tK matrix is given by:

$${}^{t}K=\left[\begin{array}{cc}\frac{\partial \text{(}{K}_{nl11}\text{(}u\text{)}+K_{nl12}\text{(}u,w\text{))}}{\partial u}& \frac{\partial \text{(}{K}_{nl11}\text{(}u\text{)}+K_{nl12}\text{(}u,w\text{))}}{\partial w}\\ \frac{\partial \text{(}{K}_{nl22}\text{(}w\text{)}+K_{nl21}\text{(}u,w\text{))}}{\partial u}& \frac{\partial \text{(}{K}_{nl22}\text{(}w\text{)}+K_{n\text{l21}}\text{(}u,w\text{))}}{\partial w}\end{array}\right]$$(A.7)

which results in:

$${}^{t}K=\left[K_l+2K_{nl}+K^{*}+K^{\text{Δ}t}\right]$$(A.8)

where:

$$\begin{array}{l}{K}^{\text{Δ}t}\left(1,1\right)=E\left(\begin{array}{l}I{w}^T{B}_{w,yy}^{nl}{}^T{B}_{u,y}+I{u}^T{B}_{u,yy}^{nl}{}^T{B}_{u,y}+I{w}^T{B}_{w,y}^{nl}{}^T{B}_{u,yy}+\dots \\ +I{u}^T{B}_{u\hspace{0.17em},y}^{nl}{}^T{B}_{u,yy}+w^T{B}_{w,y}^{nl}{}^T{B}_{u,y}+u^T{B}_{u,y}^{nl}{}^T{B}_{u,y}\end{array}\right)+\dots \\ \dots +G\left(2I{w}^T{B}_{w,y}^{nl}{}^T{B}_{u,yy}+2I{u}^T{B}_{u,y}^{nl}{}^T{B}_{u,yy}+2I{w}^T{B}_{w,yy}^{nl}{}^T{B}_{u,yy}+2I{u}^T{B}_{u,yy}^{nl}{}^T{B}_{u,yy}\right)\\ K^{\text{Δ}t}\left(2,2\right)=E\left(\begin{array}{l}I{w}^T{B}_{w,yy}^{nl}{}^T{B}_{w,y}+I{u}^T{B}_{u,yy}^{nl}{}^T{B}_{w,y}+I{w}^T{B}_{w,y}^{nl}{}^T{B}_{w,yy}+\dots \\ +I{u}^T{B}_{u\hspace{0.17em},y}^{nl}{}^T{B}_{w,yy}+w^T{B}_{w,y}^{nl}{}^T{B}_{w,y}+u^T{B}_{u,y}^{nl}{}^T{B}_{w,y}\end{array}\right)+\dots \\ \dots +G\left(2I{w}^T{B}_{w,y}^{nl}{}^T{B}_{w,yy}+2I{u}^T{B}_{u,y}^{nl}{}^T{B}_{w,yy}+2I{w}^T{B}_{w,yy}^{nl}{}^T{B}_{w,yy}+2I{u}^T{B}_{u,yy}^{nl}{}^T{B}_{w,yy}\right)\\ K^{\text{Δ}t}\left(1,2\right)=K^{\text{Δ}t}\left(2,1\right)=0\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{B}_u^{nl}=u^T{N}_u{}^T{N}_u;\hspace{1em}{B}_w^{nl}=w^T{N}_w{}^T{N}_w;\hspace{1em}{B}_u=N_u{}^T{N}_u;\hspace{1em}{B}_w=N_w{}^T{N}_w\end{array}$$

Appendix B Matrices and vectors of the equation of motion

The following matrices and vectors refer to Equation (20)

$${\text{M}}_D=\left[\begin{array}{cccc}{m}_D& 0& 0& 0\\ 0& m_D& 0& 0\\ 0& 0& I_{Dm}& 0\\ 0& 0& 0& I_{Dm}\end{array}\right]$$(B.1)

$${\text{D}}_D=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& -I_{Dy}\text{Ω}\\ 0& 0& I_{Dy}\text{Ω}& 0\end{array}\right]$$(B.2)

$${\text{D}}_D^{\text{*}}=\left[\begin{array}{cccc}0& 2m_D{\beta }_s& 0& 0\\ -2m_D{\dot{\beta }}_s& 0& 0& 0\\ 0& 0& 0& 2{\dot{\beta }}_s{I}_{Dm}-{\dot{\beta }}_s{I}_{Dy}\\ 0& 0& -2{\dot{\beta }}_s{I}_{Dm}+{\dot{\beta }}_s{I}_{Dy}& 0\end{array}\right]$$(B.3)

$$\begin{array}{l}{K}_D^{*}=[\begin{array}{ccc}-m_D\left({\dot{\gamma }}_s^2+{\dot{\beta }}_s^2\right)& m_D\left({\dot{\alpha }}_s{\dot{\gamma }}_s+{\ddot{\beta }}_s\right)& 0\\ m_D\left({\dot{\alpha }}_s{\dot{\gamma }}_s+{\ddot{\beta }}_s\right)& -m_D\left({\dot{\alpha }}_s^2+{\dot{\beta }}_s^2\right)& 0\\ 0& 0& \left(I_{Dm}-I_{Dy}\right)\left({\dot{\gamma }}_s^2-{\dot{\beta }}_s^2\right)+I_{Dy}{\dot{\beta }}_s\text{Ω}\\ 0& 0& \left({\dot{\alpha }}_s{\dot{\gamma }}_s+{\ddot{\beta }}_s\right)\left(I_{Dy}-I_{Dm}\right)\end{array}\\ \begin{array}{c}0\\ 0\\ I_{Dy}{\dot{\alpha }}_s{\dot{\gamma }}_s+\left(-{\dot{\alpha }}_s{\dot{\gamma }}_s+{\ddot{\beta }}_s\right)I_{Dm}\end{array}]\end{array}$$(B.4)

$$F_D^{*}=\left[\begin{array}{c}{m}_D\left[\ddot{X}-2{\dot{\gamma }}_s\dot{Y}+2{\dot{\beta }}_s\dot{Z}-\left({\dot{\gamma }}_s^2+{\dot{\beta }}_s^2\right)X+\left({\dot{\alpha }}_s{\dot{\beta }}_s-{\ddot{\gamma }}_s\right)\left(Y+y\right)+\left({\ddot{\beta }}_s+{\dot{\gamma }}_s{\dot{\alpha }}_s\right)Z\right]\\ m_D\left[\ddot{Z}-2{\dot{\beta }}_s\dot{X}+2{\dot{\alpha }}_s\dot{Y}-\left({\ddot{\beta }}_s-{\dot{\alpha }}_s{\dot{\gamma }}_s\right)X+\left({\ddot{\alpha }}_s-{\dot{\beta }}_s{\dot{\gamma }}_s\right)\left(Y+y\right)-\left({\dot{\beta }}_s^2+{\dot{\alpha }}_s^2\right)Z\right]\\ \begin{array}{l}{I}_{Dm}\left({\ddot{\alpha }}_s-{\dot{\beta }}_s{\dot{\gamma }}_s\right)-I_{Dm}{\dot{\gamma }}_s\left({\dot{\beta }}_s+\text{Ω}\right)\\ I_{Dm}\left({\ddot{\gamma }}_s-{\dot{\alpha }}_s{\dot{\beta }}_s\right)-I_{Dm}{\dot{\alpha }}_s\left({\dot{\beta }}_s+\text{Ω}\right)\end{array}\end{array}\right]$$(B.5)

where $I_{D_m}=\left(I_{D_x}+I_{D_z}\right)/2$

$${\text{M}}_E=\rho S{\text{M}}_1+\rho I_m{\text{M}}_2$$(B.6)

$$D_E=2\rho \text{(}{\dot{\beta }}_{s}S{D}_1-I_m\text{Ω}{D}_2\text{)}$$(B.7)

$$K_E^{\text{*}}=\rho S\left(K_1\text{+}{\ddot{\beta }}_s{D}_1+{\dot{\gamma }}_s{\dot{\alpha }}_s{K}_2\right)+\rho I_m\left[K_3-\left({\ddot{\beta }}_s+{\dot{\gamma }}_s{\dot{\alpha }}_s\right)M_4+2{\dot{\beta }}_s\left({\dot{\beta }}_s+\text{Ω}\right)M_2\right]$$(B.8)

$$\begin{array}{l}{F}_E^{\text{*}}=\rho S\left[\ddot{X}+2{\dot{\beta }}_s\dot{Z}-2{\dot{\gamma }}_s\dot{Y}+\left({\ddot{\beta }}_s+{\dot{\gamma }}_s{\dot{\alpha }}_s\right)Z+\left({\dot{\alpha }}_s{\dot{\beta }}_s-{\ddot{\gamma }}_s\right)Y-\left({\ddot{\gamma }}_s^2\text{+}{\dot{\beta }}_s^2\right)X\right]V_1\dots \\ \hspace{1em}\text{+}\left[\ddot{Z}+2{\dot{\alpha }}_s\dot{Y}-2{\dot{\beta }}_s\dot{X}+\left({\dot{\alpha }}_s{\dot{\gamma }}_s-{\ddot{\beta }}_s\right)X+\left({\ddot{\alpha }}_s+{\dot{\gamma }}_s{\dot{\beta }}_s\right)Y-\left({\dot{\alpha }}_s^2+{\dot{\beta }}_s^2\right)Z\right]V_2\dots \\ \hspace{1em}+\left({\dot{\alpha }}_s{\ddot{\beta }}_s-\ddot{\gamma }\right)V_3+\left({\ddot{\alpha }}_s+{\dot{\gamma }}_s{\dot{\beta }}_s\right)V_4\dots \\ \hspace{1em}+\rho I_m\left\{\left({\dot{\alpha }}_s{\dot{\beta }}_s-{\ddot{\gamma }}_s\right)V_5+\left({\ddot{\alpha }}_s+{\dot{\gamma }}_s{\dot{\beta }}_s\right)V_6-2\left({\dot{\beta }}_s+\text{Ω}\right)\left({\dot{\alpha }}_s{V}_5+{\dot{\gamma }}_s{V}_6\right)\right\}\end{array}$$(B.9)

$$M_1=\frac{L}{420}\left[\begin{array}{cccccccc}156& 0& 0& -22L& 54& 0& 0& 13L\\ & 156& 22L& 0& 0& 54& -13L& 0\\ & & 4L^2& 0& 0& 13L& -3L^2& 0\\ & & & 4L^2& -13L^2& 0& 0& -3L^2\\ & & & & 156& 0& 0& 22L\\ & SIM.& & & & 156& -22L& 0\\ & & & & & & 4L^2& 0\\ & & & & & & & 4L^2\end{array}\right]$$(B.10)

$$M_2=\frac{L}{420}\left[\begin{array}{cccccccc}36& 0& 0& -3L& -36& 0& 0& -3L\\ & 36& 3L& 0& 0& -36& 3L& 0\\ & & 4L^2& 0& 0& 13L& -L^2& 0\\ & & & 4L^2& 3L& 0& 0& -L^2\\ & & & & 36& 0& 0& 3L\\ & SIM.& & & & 36& -3L& 0\\ & & & & & & 4L^2& 0\\ & & & & & & & 4L^2\end{array}\right]$$(B.11)

$$M_3=\frac{1}{30L}\left[\begin{array}{cccccccc}0& 36& 3L& 0& 0& -36& 3L& 0\\ & 0& 0& -3L& -36& 0& 0& -3L\\ & & 0& -4L^2& -3L& 0& 0& L^2\\ & & & 0& 0& 3L& L^2& 0\\ & & & & 0& 36& -3L& 0\\ & SIM.& & & & 0& 0& 3L\\ & & & & & & 0& -4L^2\\ & & & & & & & 0\end{array}\right]$$(B.12)