An outliers-free isogeometric modeling method of rotating disk-shaft systems under elastic boundary conditions

Abstract

An outliers-free isogeometric modeling method for rotating disk-shaft systems is developed. The Timoshenko beam theory and artificial spring technique are employed for the rotating shaft and elastic boundary conditions. The nonlinear parameterization method is employed for the removal of outliers and three different nonlinear mappings are developed for the discussion of the accuracy of low modes. The energy coupling method between disks and shaft under nonlinear mapping is performed by using the Newton Raphson method. The results show that the isoparametric mapping has better performance in the accuracy of low modes than other nonlinear mapping and the outliers can also be removed, besides, the present method has good convergence rate for different boundary conditions. The accuracy of the proposed method shows good consistency with the Finite Element Method. The time cost of modeling is reduced by 71.4% compared to the traditional rotor model for a multiple disks rotor system, which indicates that the present approach has potential to provide more efficient optimization models of disk-shaft systems. The proposed method can provide a new modeling framework and can be easily extended to the prediction and optimization of vibration characteristics of complex rotor systems with multiple disks and supports.

Graphical abstract

Keywords

isogeometric analysis / outliers-free / non-linear mapping / disk-shaft coupling / rotor systems

1 Introduction

Rotating disk-shaft system is usually supported by elastic structures such as bearings. All of these are generally called rotor system, which are widely used in rotating machinery. With the increase in requirements for the vibration performance of modern rotating machinery, the structural complexity of rotor systems has significantly increased. The modeling method has become the most critical issue in predicting and optimizing the vibration characteristics of complex rotor systems.

Currently, the main modeling method of rotor systems is Finite Element Method (FEM). FEM was used for the first time for rotor systems by Ruhl and Booker [1]. Subsequently, the finite element model of rotor systems was improved by considering gyroscopic and shear effects [2,3]. Based on this model, there has been extensive research regarding the prediction and optimization of vibration characteristics of complex rotor systems [4–11]. Although there are many advantages for FEM in dealing with complex rotor systems [1], significant human intervention is required to rebuild the discretization model due to the difference of basis function between the finite element model and Computer Aided Design (CAD) model [12]. Therefore, a gap exists between CAD and Computer Aided Engineering (CAE) [13].

One way to eliminate the gap is Isogeometric Analysis (IGA) [14]. IGA is a discretization method developed in 2005. Its basis function, Non-Uniform Rational B-Splines (NURBS) [15], is the same with CAD systems. The discretization model can be obtained directly by the NURBS basis function and control points in a CAD model. Consequently, the gap can be eliminated and the modeling steps can be simplified. This method can produce more accurate results than corresponding finite elements in structural vibration such as thin beams [16]. IGA can also provide highly accurate solutions for spinning composite shafts [17], but the outlier modes have not been discussed in this work.

The outlier modes often appear when using IGA in structural vibration. Many methods [16,18–20] has been proposed to eliminate the outliers. Most of these methods can be roughly divided into two types: 1) changing the basis functions, such as reduced spline spaces method [19]; 2) changing the parameterization, such as the nonlinear parameterization method [16]. The former creates a new basis in reduced spline spaces by imposing additional boundary conditions, and the latter is more direct by using uniformly spaced control points to change the Jacobian. The nonlinear parameterization method will be used to remove the outliers in our work because the uniformly spaced control points can be easily achieved in the rotor system without additional boundary conditions. Besides, the lower accuracy of low modes [20] induced by the nonlinear parameterization method is also discussed in our work and a special mapping is found to be better than other nonlinear mappings for the low modes.

In view of above discussion, an outliers-free isogeometric modeling method with high accuracy in low modes for rotating disk-shaft systems is developed. The paper is structured as follows. In Section 2, the Timoshenko beam theory and artificial spring technique are employed for the rotating shaft and elastic boundary conditions, and Hamilton’s principle is used to derive the motion equation. In Section 3, three different nonlinear mappings based on the nonlinear parameterization method are provided for the discussion of the accuracy of low modes, and the approach of the energy coupling between disks and shaft under a nonlinear mapping is proposed by using Newton Raphson method. In Section 4, the accuracy of low modes, the removal of outliers and the convergence rate of a Timoshenko beam with different boundary conditions are discussed; moreover, the accuracy of the isogeometric rotor system with different parameters is examined and the further validation in multiple disks rotor system is presented, the characteristics and advantages of the present method are also discussed. Finally, in Section 5, the conclusions and future research work are given.

2 Theory of the isogeometric rotor system

In this section, the potential energy and kinetic energy of rotor system are derived based on Timoshenko beam theory and artificial spring technique. The Hamilton’s principle is employed to obtain the discretized form of the equation of motions.

2.1 The structure of rotor system

Fig.1 shows the rotor system studied here, where an annular disk with mass m_d, the polar moment of inertia, J_p, and the diametral moment of inertia, J_d, is mounted on a shaft with radius r and length L. The distance from the shaft left end to the disk is a. The disk-shaft system is supported by two elastic supports using artificial spring technique [21] and each support has principal stiffnesses, k_xx and k_yy, and Cross-coupled stiffnesses, k_xy and k_yx. The inertial coordinate system xyz is defined and the x-axis aligned along the undeformed shaft centerline direction.

Fig.1 Rotor system and coordinate system.

Full size|PPT slide

2.2 Energy equations of the rotor system

The normal strain ε_zzx, ε_zzy, and shear strain γ_zx and γ_zy of a shaft are obtained by Timoshenko beam theory,

(1)$\begin{array}{l}{\text{ε}}_{\mathit{\text{zz}}x}=u{\text{θ}}_{y,\mathit{\text{z}}},{\text{ε}}_{\mathit{\text{zz}}y}=v{\text{θ}}_{x,\mathit{\text{z}}},\\ {\text{γ}}_{\mathit{\text{z}}x}={\text{θ}}_y-{\mathit{u}}_{,\mathit{\text{z}}},{\text{γ}}_{\mathit{\text{z}}y}={\text{θ}}_x-{\mathit{v}}_{,\mathit{\text{z}}},\end{array}$

where u and v denote the translations of the shaft with respect to x-and y-axes, respectively. θ_x and θ_y are the rotations of the shaft with respect to x-and y-axes, respectively. The potential energy of the shaft induced by the bending deformation is given by

(2)$\begin{array}{rl}\delta U_{\text{s}}=& \underset{V}{\int }(\delta {\text{ε}}_{\mathit{\text{zz}}x}^{\text{T}}E{\text{ε}}_{\mathit{\text{zz}}x}+\delta {\text{ε}}_{\mathit{\text{zz}}y}^{\text{T}}E{\text{ε}}_{\mathit{\text{zz}}y}\\ & +\delta {\text{γ}}_{\mathit{\text{z}}x}^{\text{T}}\kappa G\delta {\text{γ}}_{\mathit{\text{z}}x}+\delta {\text{γ}}_{\mathit{\text{z}}y}^{\text{T}}\kappa G\delta {\text{γ}}_{\mathit{\text{z}}y})\mathrm{d}V.\end{array}$

where E denotes Young’s modulus, κ and G are shear coefficient and shear modulus, respectively.

The motion of a material point P on the arbitrary cross section of a rotating shaft can be divided into two parts, translations and rotations (Fig.2). The position vector of point P with respect to the coordinate OXYZ induced by translations is

Fig.2 The coordinate systems for the description of the motion of a rotary shaft: (a) description of the translations; (b) description of the rotations.

Full size|PPT slide

(3)$\mathit{O}\mathit{P}=\mathit{O}\mathit{C}+\mathit{C}\mathit{P}+\left\{\begin{array}{c}u\\ v\\ w\end{array}\right\}.$

Therefore, the velocity of translations can be defined as

(4)${\mathit{v}}_{\text{ts}}=\frac{\text{d}(\mathit{O}\mathit{P})}{\text{d}t}=\left\{\begin{array}{c}\dot{u}\\ \dot{v}\\ \dot{w}\end{array}\right\},$

where w is set to 0 because we focus on the transverse motion. To describe the rotations of the point, suppose the rotations of coordinate are applied in the following order (Fig.2(b)), θ_y about y-axis obtaining new coordinate ox₁y₁z₁, θ_x about x₁-axis obtaining new coordinate ox₂y₂z₂, and then θ_z = Ωt about z₂-axis obtaining coordinate ox₃y₃z₃, where Ω is the rotation speed of the shaft. The instantaneous velocity of the point P is expressed as,

(5)${\mathit{v}}_{\text{rs}}={\text{ω}}_{\text{s3}}\times \mathit{r},$

where r = {x₃, y₃, 0}^T denotes the position vector of point P related to the coordinate ox₃y₃z₃, and the angular velocity can be expressed by

(6)${\text{ω}}_{\text{s3}}=\left\{\begin{array}{c}{\dot{\theta }}_{x_3}\\ {\dot{\theta }}_{y_3}\\ {\dot{\theta }}_{{\mathit{\text{z}}}_3}\end{array}\right\}={\text{T}}_{\text{r}}\left\{\begin{array}{c}{\dot{\theta }}_x\\ {\dot{\theta }}_y\\ {\dot{\theta }}_{\mathit{\text{z}}}\end{array}\right\},$

where T_r is the transformation matrix defined by

(7)$T_{\text{r}}=\left[\begin{array}{ccc}\cos {\theta }_{\mathit{\text{z}}}& \sin {\theta }_{\mathit{\text{z}}}\cos {\theta }_x& 0\\ -\sin {\theta }_{\mathit{\text{z}}}& \cos {\theta }_{\mathit{\text{z}}}\cos {\theta }_x& 0\\ 0& -\sin {\theta }_x& 1\end{array}\right].$

The variation of ω_s3 can be derived as

(8)$\begin{array}{rl}\delta {\omega }_{\text{s3}}=& [\left(-\sin {\theta }_{\mathit{\text{z}}}{\dot{\theta }}_x+\cos {\theta }_{\mathit{\text{z}}}\cos {\theta }_x{\dot{\theta }}_y\right)\delta {\theta }_{\mathit{\text{z}}}+\cos {\theta }_{\mathit{\text{z}}}\delta {\dot{\theta }}_x\\ & -\sin {\theta }_{\mathit{\text{z}}}\sin {\theta }_x{\dot{\theta }}_y\delta {\theta }_x+\sin {\theta }_{\mathit{\text{z}}}\cos {\theta }_x\delta {\dot{\theta }}_y-\\ & \left(\cos {\theta }_{\mathit{\text{z}}}{\dot{\theta }}_x-\sin {\theta }_{\mathit{\text{z}}}\cos {\theta }_x{\dot{\theta }}_y\right)\delta {\theta }_{\mathit{\text{z}}}-\sin {\theta }_{\mathit{\text{z}}}\delta {\dot{\theta }}_x\\ & -\cos {\theta }_{\mathit{\text{z}}}\sin {\theta }_x{\dot{\theta }}_y\delta {\theta }_x+\cos {\theta }_{\mathit{\text{z}}}\cos {\theta }_x\delta {\dot{\theta }}_y\\ & -\cos {\theta }_x{\dot{\theta }}_y\delta {\theta }_x-\sin {\theta }_x\delta {\dot{\theta }}_y+\delta {\dot{\theta }}_{\mathit{\text{z}}}].\end{array}$

The kinetic energy of the rotating shaft can be derived by the formula:

(9)$\delta T_{\text{s}}=\rho \int (\delta {\mathit{v}}_{\text{ts}}^{\text{T}}{\mathit{v}}_{\text{ts}}^{}+\delta {\mathit{v}}_{\text{rs}}^{\text{T}}{\mathit{v}}_{\text{rs}})\text{d}V.$

The kinetic energy of the rotating disk, δT_d, can be obtained in the same way. For the convenience of discussion, the moment of inertia ratio δ of disk is defined as

(10)$\delta =\frac{J_{\text{p}}}{J_{\text{d}}}.$

In general elastic supports, the principal stiffness is enough to obtain accuracy results of rotor system. But in some special elastic supports such as gas/liquid seals, hydrodynamic bearings, the Cross-coupled stiffness need to be considered [22]. Therefore, the elastic boundary condition at arbitrary position z = z_b can be obtained by using artificial spring technology,

(11)$\delta U_{\text{b}}={{\left\{\begin{array}{c}\delta u\\ \delta v\end{array}\right\}}^{\text{T}}{\mathit{K}}_{\text{b}}\left\{\begin{array}{c}u\\ v\end{array}\right\}|}_{\mathit{\text{z}}={\mathit{\text{z}}}_{\text{b}}},$

where

(12)${\mathit{K}}_{\text{b}}=\left[\begin{array}{cc}-k_{xx}& -k_{xy}\\ -k_{yx}& -k_{yy}\end{array}\right].$

For the convenience of discussion, we defined

(13)$k_{yy}=\alpha k_{xx},k_{xy}={\beta }_1{k}_{xx},k_{yx}={\beta }_2{k}_{xx},$

where α, β₁, β₂ are dimensionless coefficients.

2.3 Free vibration analysis

The weak form of the equations of motion can be derived by using Hamilton’s principle,

(14)${\int }_{t_2}^{t_2}(\delta T-\delta U)\text{d}t=0,$

where

(15)$\begin{array}{c}\delta T=\delta T_{\text{s}}+\delta T_{\text{d}},\\ \delta U=\delta U_{\text{s}}+\delta U_{\text{b}}.\end{array}$

By using NURBS-based approximation function, the IGA discretized form of the equation of motions for the vibration of rotor system can be represented as

(16)$({\mathit{M}}_{\text{s}}+{\mathit{M}}_{\text{d}})\ddot{\mathit{q}}+({\mathit{G}}_{\text{s}}+{\mathit{G}}_{\text{d}})\dot{\mathit{q}}+({\mathit{K}}_{\text{s}}+{\mathit{K}}_{\text{b}})\mathit{q}=0,$

where M, G, K are mass, gyroscopic, and stiffness matrices, respectively, the subscripts s, d, b denote shaft, disks, and supports, respectively. To solve Eq. (16) using the state space method, the above equation can be rearranged as the following form,

(17)$\begin{array}{rl}& \left[\begin{array}{cc}{\mathit{G}}_{\text{s}}+{\mathit{G}}_{\text{d}}& {\mathit{M}}_{\text{s}}+{\mathit{M}}_{\text{d}}\\ {\mathit{M}}_{\text{s}}+{\mathit{M}}_{\text{d}}& 0\end{array}\right]\frac{\text{d}}{\text{d}t}\left\{\begin{array}{c}\mathit{q}\\ \dot{\mathit{q}}\end{array}\right\}\\ +& \left[\begin{array}{cc}{\mathit{K}}_{\text{s}}+{\mathit{K}}_{\text{b}}& 0\\ 0& -({\mathit{M}}_{\text{s}}+{\mathit{M}}_{\text{d}})\end{array}\right]\left\{\begin{array}{c}\mathit{q}\\ \dot{\mathit{q}}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\end{array}\right\},\end{array}$

which can be rewritten as

(18)$\mathit{A}\dot{\mathit{x}}+\mathit{B}\mathit{x}=0,$

where

(19)$\begin{array}{c}\mathit{x}=\left\{\begin{array}{c}\mathit{q}\\ \dot{\mathit{q}}\end{array}\right\},\mathit{A}=\left[\begin{array}{cc}{\mathit{G}}_{\text{s}}+{\mathit{G}}_{\text{d}}& {\mathit{M}}_{\text{s}}+{\mathit{M}}_{\text{d}}\\ {\mathit{M}}_{\text{s}}+{\mathit{M}}_{\text{d}}& 0\end{array}\right],\\ \mathit{B}=\left[\begin{array}{cc}{\mathit{K}}_{\text{s}}+{\mathit{K}}_{\text{b}}& 0\\ 0& -({\mathit{M}}_{\text{s}}+{\mathit{M}}_{\text{d}})\end{array}\right].\end{array}$

Suppose x = x₀e^λt, by substituting this equation into Eq. (18), the solution to the eigenvalue problem can be obtained by

(20)$(\lambda \mathit{A}+\mathit{B}){\mathit{x}}_0=0,$

where λ is a complex number which includes natural frequencies and x₀ is the vector of the degrees of freedom.

3 Removal of outliers and energy coupling method

In this section, three different nonlinear mappings and the energy coupling method between disk and shaft are presented based on the nonlinear parameterization method.

3.1 Change of parameterization

The nonlinear parameterization method is proposed in Ref. [16] to remove the outliers. This method can be directly implemented in the rotor system without a change of basis. The uniformly spaced control points need to be defined first. For a rotor system with length L, the control points can be defined as

(21)${\mathit{\text{z}}}_{\text{CP}}=\left\{0,\frac{L}{n_{\text{el}}+p_{\text{z}}},\frac{2L}{n_{\text{el}}+p_{\text{z}}},\cdots ,\frac{(n_{\text{el}}+p_{\text{z}}-1)L}{n_{\text{el}}+p_{\text{z}}},L\right\},$

where n_el denotes the number of element and p_z denotes the degree of the mapping function from physical domain to parameter domain, the mapping function is defined as following

(22)$\mathit{\text{z}}(\xi )={\mathit{R}}_{p_{\text{z}}}{\mathit{\text{z}}}_{\text{CP}},$

where R is the NURBS basis.

However, the worse accuracy of low modes limited the utilization of this method, which hasn’t been further discussed. In this paper, three different nonlinear mappings are presented for the better observation of the accuracy of low modes:

1) subparametric mapping, 1 < p_z < p_u, where p_z = 1 is ignored because we focus on nonlinear mapping;

2) isoparametric mapping, p_z = p_u;

3) superpaprametric mapping, p_z > p_u, where p_u denotes the degree of NURBS basis of displacement (u,v) and rotation (θ,φ). The accuracy of low frequencies will be further discussed in Section 4.

3.2 Energy coupling method

For the energy coupling between disk and shaft in IGA with the nonlinear parameterization method, the following steps should be taken into account.

1) Determination of parametric coordinate of coupling point.

To get the parametric coordinate ξ_d of coupling point, a nonlinear equation about ξ_d based the mapping relationship is defined as

(23)$f(\xi )=\left|\mathit{\text{z}}(\xi )-a\right|,$

whose minimum condition with respect to the parametric coordinate leads to

(24)$f(\xi )=0.$

Consider the Taylor series of Eq. (24) around ξ^(k) where k is an iteration counter.

(25)$f({\xi }^{(k+1)})=f({\xi }^{(k)})+f^{\prime }({\xi }^{(k)})({\xi }^{(k+1)}-{\xi }^{(k)})+o({\xi }^{(k+1)}).$

A Newton iterative process is invoked to find the parametric coordinate.

(26)${\xi }^{(k+1)}={\xi }^{(k)}-\frac{f({\xi }^{(k)})}{f^{{}^{\prime }}({\xi }^{(k)})}.$

Using point coincidence convergence criteria

(27)$f({\xi }^{(k+1)})⩽\epsilon .$

The converged solution determines the parametric coordinate ξ_d

(28)${\xi }_{\text{d}}={\xi }^{(k+1)}.$

2) Establishment of energy coupling functions.

The field functions of coupling point ξ_d are defined by

(29)$\{\begin{array}{c}{u}_{\text{d}}=\sum _{i=1}^{p+1}{R}_{i,p}^{(j_{\text{de}})}({\xi }_{\text{d}})u_i^{(j_{\text{de}})},\\ v_{\text{d}}=\sum _{i=1}^{p+1}{R}_{i,p}^{(j_{\text{de}})}({\xi }_{\text{d}})v_i^{(j_{\text{de}})},\\ {\phi }_{\text{d}}=\sum _{i=1}^{p+1}{R}_{i,p}^{(j_{\text{de}})}({\xi }_{\text{d}}){\phi }_i^{(j_{\text{de}})},\\ {\theta }_{\text{d}}=\sum _{i=1}^{p+1}{R}_{i,p}^{(j_{\text{de}})}({\xi }_{\text{d}}){\theta }_i^{(j_{\text{de}})},\end{array}$

where u_d, v_d, θ_d, and φ_d are the translations and rotations of the coupling point, j_de denotes that the disk is at the j_deth beam element. The matrix of disk can be obtained by substituting Eq. (29) into δT_d.

Those steps are applicable for general disk-shaft coupling but explained for support-shaft coupling without loss of generality.

4 Results and discussion

In this section, the low frequency error and L² error of mode shape of Timoshenko beam with different boundary conditions and different degree of basis are presented. Besides, the whole frequency spectrum and convergency rate of Timoshenko beam with different boundary conditions are also presented. Furthermore, the accuracy of rotor systems with different parameters and the characteristic and advantage of isogeometric rotor system are discussed by comparing with FEM.

4.1 Accuracy of low modes, outliers and convergence rate

Fig.3 and Fig.4 illustrate the low frequency error and L² error of mode shape of Timoshenko beam with simply supported and free boundary conditions. The material and dimension parameters are listed in Tab.1. Every model has a fixed degree of freedom, N = 200. The degree p_z ranges from 2 to 4 for p_u = 3, 3 to 5 for p_u = 4, and 4 to 6 for p_u = 5, so that the three mappings can all be included. The error of simply supported and free boundary condition is calculated by using the analytical solution from Ref. [23] and a very fine mesh with standard B-splines, respectively. It can be observed that the accuracy of the first ten frequencies and mode shapes obtained by isoparametric mapping is always better than subparametric and superparametric mapping in all p_u and different boundary conditions. The reason is that the Jacobians of the three mappings are different, which leads to the different performance in low modes. On the other hand, the outliers can also be removed by isoparametric mapping, see Fig.5 and Fig.6. More detail of the removal of outliers using nonlinear parameterization method can be found in Ref. [16]. The above results show that the isoparametric mapping has better performance in the accuracy of low modes by comparing with subparametric and superparametric mapping, and the outliers can also be improved.

Tab.1 the parameters of rotor system

E (GPa)	µ	ρ (kg/m3)	L (m)	r (m)	md (kg)	Jd (kg∙m2)	kxx (N/m)
210	0.3	7850	1	0.02	40	0.1	1.00 × 107

Fig.3 The first ten Frequency error (left) and L2 error of the mode shapes (right) of free vibration of a simply supported Timoshenko beam: (a) pu = 3; (b) pu = 3; (c) pu = 4; (d) pu = 4; (e) pu = 5; (f) pu = 5.

Full size|PPT slide

Fig.4 The first ten Frequency error (left) and L2 error of the mode shapes (right) of free vibration of a free-free Timoshenko beam: (a) pu = 3; (b) pu = 3; (c) pu = 4; (d) pu = 4; (e) pu = 5; (f) pu = 5.

Full size|PPT slide

Fig.5 The whole frequency spectrum of free vibration of a simply supported Timoshenko beam.

Full size|PPT slide

Fig.6 The whole frequency spectrum of free vibration of a free-free Timoshenko beam.

Full size|PPT slide

Fig.7–Fig.9 illustrate the convergence of the first four forward frequencies of a simply supported spinning Timoshenko shaft at different rotation speeds, based on the reference solution that can be found in Ref. [24]. The results show that the model has good convergence at different rotation speeds.

Fig.7 The convergence rate of the first four forward frequencies of a simply supported spinning shaft, 0 r/min, cubic basis.

Full size|PPT slide

Fig.8 The convergence rate of the first four forward frequencies of a simply supported spinning shaft, 500 r/min, cubic basis.

Full size|PPT slide

Fig.9 The convergence rate of the first four forward frequencies of a simply supported spinning shaft, 1000 r/min, cubic basis.

Full size|PPT slide

4.2 Accuracy of isogeometric rotor system in different parameters

A rotor system includes shaft, disks and elastic supports and so on, whose complex structures make it hard to obtain analytical solutions. Therefore, the accuracy of the isogeometric rotor system is compared with FEM in different parameters, including different boundary conditions (cantilever, simply support), different parameters of supports (isotropic, anisotropic, cross-coupled stiffnesses) and offset disk, and the further validation in multiple disks rotor system is also presented.

To verify the accuracy of the present method for different parameters, five rotor systems based on Fig.1 with different boundary conditions, different parameters of supports and offset disk are considered.

1) Case 1: elastic boundary conditions are at both shaft ends; without cross-coupled stiffness; isotropic principal stiffness is considered. δ = 2, a = 0.5, β₁ = 0, β₂ = 0, α = 1.

2) Case 2: elastic boundary conditions are at both shaft ends; without cross-coupled stiffness; anisotropic principal stiffness is considered. δ = 2, a = 0.5, β₁ = 0, β₂ = 0, α = 0.5.

3) Case 3: elastic boundary conditions are at both shaft ends; anisotropic cross-coupled stiffness and principal stiffness are considered. δ = 2, a = 0.5, β₁ = 0.5, β₂ = 0.5, α = 0.5.

4) Case 4: elastic boundary conditions are at the left end of shaft and the right end of shaft is free; without cross-coupled stiffness; isotropic principal stiffness is considered; the disk is located at the right end of shaft. δ = 2, a = 1, β₁ = 0, β₂ = 0, α = 1.

5) Case 5: elastic boundary conditions are at both shaft ends; offset disk; without cross-coupled stiffness; isotropic principal stiffness is considered. δ = 2, a = 0.3, β₁ = 0, β₂ = 0, α = 1.

The first eight natural frequencies of the five cases are compared in Tab.2–Tab.6, both isogeometric model and finite element (FE) model have 252 degrees of freedom and the degree of basis is 3. It is found that the results are consistent with FEM, and maximum error is found to be 1.34%, which indicates the present method can handle different rotor systems and achieve accurate results. The errors might be caused by the difference of the coupling form of the disk. In FEM, the disk is coupled at one node, and the kinetic energy is centralized at the node. However, the actual kinetic energy of the disk is distributed at an interval whose width might be the thickness of disk. The distribution of the kinetic energy might be considered in IGA because the disk is coupled at one element. it might provide a possibility to simulated the actual kinetic energy distribution of the disk by using IGA. On the other hand, it is also found that the results of Case 4 are more consistent with FEM, the reason is the coupling form of the disk at both shaft ends is the same with FEM due to the endpoint interpolation property of NURBS basis (Ref. [25], Subsection 3.2).

Tab.2 Accuracy of the natural frequencies (Hz) of Case 1, δ = 2, a = 0.5, β1 = 0, β2 = 0, α = 1

Mode	0 r/min				500 r/min				1000 r/min
	FEM	Present	Error		FEM	Present	Error		FEM	Present	Error
1	6.962	6.960	0.03%		6.962	6.960	0.03%		6.962	6.960	0.03%
2	6.962	6.960	0.03%		6.962	6.960	0.03%		6.962	6.960	0.03%
3	66.522	66.953	0.65%		59.582	60.033	0.76%		53.421	53.883	0.86%
4	66.522	66.953	0.65%		74.241	74.641	0.54%		82.700	83.059	0.43%
5	249.794	249.806	0.00%		249.786	249.797	0.00%		249.777	249.789	0.00%
6	249.794	249.806	0.00%		249.803	249.815	0.00%		249.812	249.823	0.00%
7	261.267	264.458	1.22%		260.450	263.607	1.21%		259.724	262.848	1.20%
8	261.267	264.458	1.22%		262.192	265.417	1.23%		263.247	266.506	1.24%

Tab.3 Accuracy of the natural frequencies (Hz) of Case 2, δ = 2, a = 0.5, β1 = 0, β2 = 0, α = 0.5

Mode	0 r/min				500 r/min				1000 r/min
	FEM	Present	Error		FEM	Present	Error		FEM	Present	Error
1	6.948	6.946	0.03%		6.948	6.946	0.03%		6.948	6.946	0.03%
2	6.962	6.960	0.03%		6.962	6.960	0.03%		6.962	6.960	0.03%
3	66.323	66.749	0.64%		59.490	59.939	0.76%		53.336	53.796	0.86%
4	66.522	66.953	0.65%		74.132	74.529	0.54%		82.579	82.935	0.43%
5	244.902	244.912	0.00%		244.902	244.912	0.00%		244.901	244.912	0.00%
6	249.794	249.806	0.00%		249.794	249.806	0.00%		249.794	249.806	0.00%
7	256.300	259.337	1.19%		256.206	259.237	1.18%		255.958	258.972	1.18%
8	261.267	264.458	1.22%		261.471	264.668	1.22%		262.052	265.267	1.23%

Tab.4 Accuracy of the natural frequencies (Hz) of Case 3, δ = 2, a = 0.5, β1 = 0.5, β2 = 0.5, α = 0.5

Mode	0 r/min				500 r/min				1000 r/min
	FEM	Present	Error		FEM	Present	Error		FEM	Present	Error
1	6.904	6.901	0.03%		6.904	6.901	0.03%		6.904	6.901	0.03%
2	6.965	6.963	0.03%		6.965	6.963	0.03%		6.965	6.963	0.03%
3	65.677	66.087	0.62%		59.203	59.645	0.75%		53.079	53.532	0.85%
4	66.569	67.002	0.65%		73.814	74.205	0.53%		82.217	82.566	0.42%
5	229.460	229.468	0.00%		229.460	229.468	0.00%		229.460	229.468	0.00%
6	240.552	243.118	1.07%		240.577	243.143	1.07%		240.655	243.216	1.06%
7	250.950	250.962	0.00%		250.950	250.962	0.00%		250.950	250.962	0.00%
8	262.437	265.665	1.23%		262.526	265.755	1.23%		262.799	266.031	1.23%

Tab.5 Accuracy of the natural frequencies (Hz) of Case 4, δ = 2, a = 1, β1 = 0, β2 = 0, α = 1

Mode	0 r/min				500 r/min				1000 r/min
	FEM	Present	Error		FEM	Present	Error		FEM	Present	Error
1	28.228	28.230	0.01%		23.631	23.638	0.03%		19.793	19.803	0.05%
2	28.228	28.230	0.01%		33.363	33.359	0.01%		38.581	38.573	0.02%
3	77.443	77.553	0.14%		74.958	75.065	0.14%		73.141	73.245	0.14%
4	77.443	77.553	0.14%		80.890	81.000	0.14%		85.650	85.757	0.12%
5	207.993	208.262	0.13%		207.660	207.928	0.13%		207.374	207.640	0.13%
6	207.993	208.262	0.13%		208.381	208.653	0.13%		208.844	209.118	0.13%
7	421.283	421.797	0.12%		421.188	421.704	0.12%		421.102	421.615	0.12%
8	421.283	421.797	0.12%		421.381	421.897	0.12%		421.488	422.004	0.12%

Tab.6 Accuracy of the natural frequencies (Hz) of Case 5, δ = 2, a = 0.3, β1 = 0, β2 = 0, α = 1

Mode	0 r/min				500 r/min				1000 r/min
	FEM	Present	Error		FEM	Present	Error		FEM	Present	Error
1	8.204	8.200	0.05%		8.130	8.126	0.05%		8.053	8.048	0.05%
2	8.204	8.200	0.05%		8.276	8.272	0.05%		8.345	8.342	0.04%
3	69.807	70.452	0.92%		63.732	64.457	1.14%		58.148	58.929	1.34%
4	69.807	70.452	0.92%		76.270	76.810	0.71%		82.945	83.361	0.50%
5	140.645	142.100	1.03%		139.154	140.533	0.99%		137.975	139.287	0.95%
6	140.645	142.100	1.03%		142.563	144.101	1.08%		145.069	146.178	0.76%
7	410.289	412.840	0.62%		410.186	412.707	0.61%		410.102	411.759	0.40%
8	410.289	412.840	0.62%		410.392	412.984	0.63%		410.514	412.310	0.44%

To further validate the present method, the vibration characteristics of a three-disk rotor system (Fig.10) is compared with FEM, both isogeometric model and FE model have 252 degrees of freedom and the degree of basis is 3. Each disk with thickness 10 mm and diameter 199 mm is mounted on a shaft with diameter 20 mm. The material properties of the disk-shaft system are: E = 210 GPa, ρ = 7850 kg/m³, Poisson’s ratio υ = 0.3. The stiffnesses of the support at the left side of the shaft are K1_xx = K1_yy = 2 × 10⁸ N/m and the stiffnesses of the support at the right side of the shaft are K2_xx = K2_yy = 2.8 × 10⁷ N/m. As seen from Campbell diagram in Fig.11, the results obtained from the present method are consistent with FEM, which signifies that this method can accurately predict the separate phenomenon induced by gyroscopic effects. The comparison of the first two critical speeds and the corresponding mode shapes are shown in Tab.7 and Fig.12, respectively, the results show that the vibration characteristics can be obtained accurately by this method and it can be easily extended to multiple disks rotor systems.

Tab.7 The comparison of the first two critical speeds (r/min)

Mode	IGA	FEM	Difference (%)
1	5660.1	5659.5	0.01
2	13135.7	13111	0.19

Fig.10 The structure and dimensions of a three-disk rotor system.

Full size|PPT slide

Fig.11 Campbell diagram of three-disk rotor system.

Full size|PPT slide

Fig.12 The comparison of the corresponding mode shapes of the first two critical speeds: (a) 1st Mode shape; (b) 2nd Mode shape.

Full size|PPT slide

4.3 The characteristics and advantages of isogeometric rotor systems

4.3.1 Characteristics

One of the characteristics of isogeometric rotor systems is that the information of the locations for disks is not contained in shaft mesh, which is one of the differences between the present method and traditional rotor model based on standard FEM. Based on the rotor system in Fig.1, the comparison of the coarsest shaft mesh using the two methods is illustrated as shown in Fig.13. The red points and the blue point are the endpoints of elements and the blue point also represents the location of the disk. In Fig.13(a), at least two elements (three nodes) are needed in the FE shaft mesh, and the middle node must be used for the description of the location of the disk. However, only one element is needed in the shaft mesh using IGA (Fig.13(b)), and there is no need to create a node/knot to describe the location of the disk. It is a common phenomenon in isogeometric rotor systems. The reason for this phenomenon is that the discretization of shaft (h-refinement, see Subsubsection 3.2.1) was finished before using the location information of disk (determination of parametric coordinate of coupling point, see Subsubsection 3.2.2). Therefore, the information of the locations for disks doesn’t exist in isogeometric shaft mesh. To further illustrate the characteristics of isogeometric rotor systems, the general modeling processes of disk-shaft systems using standard FEM and IGA are depicted in Fig.14 and Fig.15, respectively. The left sides are the descriptions for the modeling steps and the corresponding diagrams are presented at the right sides.

Fig.13 The discretization model of shaft: (a) traditional rotor model based on FEM; (b) IGA.

Full size|PPT slide

Fig.14 The traditional disk-shaft system FE model.

Full size|PPT slide

Fig.15 The isogeometric disk-shaft system.

Full size|PPT slide

4.3.2 Advantages

The locations of disks are usually optimized in the design of rotor dynamics [26], and this optimization process can be improved based on the above characteristics. When standard FEM is employed for the optimization of the locations of disks, the traditional way including two steps [26]. 1) The regeneration of shaft mesh (Fig.16(a) and 16(b)); 2) The recalculation of the corresponding matrixes (step 1 in Fig.14). These two steps have to be repeated because the locations of disks are included in the shaft mesh and the locations of disks will change during the process of optimization. However, due to the fact that the information of the locations for disks is independent of the isogeometric shaft mesh (Fig.15), only recalculation of the coupling matrixes of disk-shaft systems needs to be repeated (step 2 in Fig.15) and the shaft mesh will not change after optimization (Fig.16(c) and Fig.16(d)).

Fig.16 The comparison of optimization model for a disk-shaft system: (a) standard FEM, before optimization; (b) standard FEM, after optimization; (c) IGA, before optimization; (d) IGA, after optimization.

Full size|PPT slide

The summarization of the modeling processes using the two methods during optimization is shown in Fig.17. The difference between these two methods is that the process for the regeneration of finite element shaft mesh and the recalculation of corresponding matrixes are replaced by the recalculation of disk-shaft coupling matrixes using IGA. Because the number of disks is generally far less than the number of shaft elements, the computational cost for the calculation of disk-shaft coupling matrixes is also less than the regeneration of shaft mesh and the calculation of the corresponding matrixes.

Fig.17 The modeling steps of rotor systems during the process of optimization: (a) standard FEM; (b) IGA.

Full size|PPT slide

A simple numerical validation using the traditional rotor model based on FEM and IGA is conducted by using the three-disk rotor system in Fig.10 with the same degrees of freedom and polynomial degree. Only the modeling process is considered because the rest optimization process is totally the same in both methods [26]. the modeling process is summarized in Fig.17. To improve the accuracy of the benchwork, the modeling processes are repeated 80 times. All processes are implemented on an Intel i7-12700KF CPU and 32 GB of memory. The comparison of time cost for modeling is shown in Fig.18. It is found that the time cost of each modeling process using IGA is less than standard FEM, and the average time cost is reduced by 71.4% compared to standard FEM, which indicates that more efficient optimization models of rotor systems can be provided by IGA. It can be conveniently extended to other situations such as the optimization of the locations of supports.

Fig.18 The comparison of time cost for the modeling of the three-disk rotor system during the process of optimization.

Full size|PPT slide

5 Conclusions

In this work we have explored IGA with outliers-free as a modeling method for disk-shaft rotor systems with different elastic boundary conditions, the energy coupling method under nonlinear mapping is proposed by using Newton Raphson method. Shear effects, gyroscopic effects, offset disk, multiple disks, different elastic boundary conditions and the different parameters of supports including isotropic/anisotropic and cross-coupled stiffnesses are considered. For the results, first, three different mappings based on nonlinear parameterization method are presented for the detailed discussion of the frequency error and mode shape L² error of low modes of Timoshenko beam with simply supported and free boundary conditions and different degree of basis. The removal of outliers and convergence rate of first four frequencies are also presented. Secondly, the accuracy of rotor systems with different parameters including different boundary conditions (cantilever, simply support), different parameters of supports (isotropic, anisotropic, cross-coupled stiffnesses), offset disk and multiple disks is presented. Finally, the characteristics and advantages of isogeometric rotor systems are also discussed. The conclusions are listed as follows.

The nonlinear parameterization method using isoparametric mapping has better performance in the accuracy of low modes by comparing with subparametric and superparametric mapping in Timoshenko model with different boundary conditions and different degree of basis funcion, and the outliers can also be improved. The present method can achieve high convergence rate around 7 in different rotation speeds.

The present approach can accurately predict the vibration characteristics of disk-shaft rotor systems in different parameters. In addition, it can be observed that the results of cantilever rotor are more in line with FEM, the reason is the coupling form of the disk at both shaft ends is the same with FEM due to the endpoint interpolation property of NURBS basis.

In the process of optimization, the regeneration of finite element shaft mesh and the recalculation of corresponding matrixes are replaced by the recalculation of disk-shaft coupling matrixes using IGA. Due to the number of disk elements being generally far less than the number of shaft elements, the computational cost for the calculation of disk-shaft coupling matrixes is also less than the regeneration of shaft mesh and the calculation of corresponding matrixes. This is validated in a three-disk rotor system and the results show that the average time cost of modeling using the developed method is reduced by 71.4% compared to traditional rotor model based on FEM. These results indicate that more efficient optimization models of rotor systems might be provided by the present method.

The proposed approach can be extended to the prediction and optimization of vibration characteristics of complex rotor systems with multiple disks and supports without difficulty, which is important to the design of complex rotor systems. However, it is worth nothing that only linear stiffnesses and kinetic energies of disks are considered by neglecting nonlinearity effects and elastic deformations of disks. In the future work, we will try to introduce the strain energy of disks and nonlinear stiffness and damping of the supports and coupling interfaces in rotor systems.

References

[1] Ruhl R L, Booker J F. A finite element model for distributed parameter turborotor systems. Journal of Engineering for Industry, 1972, 94(1): 126–132

[2] Nelson H D, McVaugh J M. The dynamics of rotor-bearing systems using finite elements. Journal of Engineering for Industry, 1976, 98(2): 593–600

[3] Nelson H D. A finite rotating shaft element using Timoshenko beam theory. Journal of Mechanical Design, 1980, 102(4): 793–803

[4] Lee A S, Kim B O, Kim Y C. A finite element transient response analysis method of a rotor-bearing system to base shock excitations using the state-space Newmark scheme and comparisons with experiments. Journal of Sound and Vibration, 2006, 297(3–5): 595–615

[5] Gayen D, Chakraborty D, Tiwari R. Whirl frequencies and critical speeds of a rotor-bearing system with a cracked functionally graded shaft––Finite element analysis. European Journal of Mechanics. A, Solids, 2017, 61: 47–58

[6] Gayen D, Tiwari R, Chakraborty D. Finite element based stability analysis of a rotor-bearing system having a functionally graded shaft with transverse breathing cracks. International Journal of Mechanical Sciences, 2019, 157: 403–414

[7] Kulesza Z, Sawicki J T. Rigid finite element model of a cracked rotor. Journal of Sound and Vibration, 2012, 331(18): 4145–4169

[8] Bhowmick A, Ganguly S, Neogy S, Nandi A. A shaft finite element for analysis of viscoelastic tapered Timoshenko rotors. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2020, 42(1): 48

[9] Belhocine A, Afzal A. Computational finite element analysis of brake disc rotors employing different materials. Australian Journal of Mechanical Engineering, 2022, 20(3): 637–650

[10] Prabith K, Krishna I R P. The numerical modeling of rotor–stator rubbing in rotating machinery: A comprehensive review. Nonlinear Dynamics, 2020, 101(2): 1317–1363

[11] Kumar P, Tiwari R. Finite element modelling, analysis and identification using novel trial misalignment approach in an unbalanced and misaligned flexible rotor system levitated by active magnetic bearings. Mechanical Systems and Signal Processing, 2021, 152: 107454

[12] RogersD F. An Introduction to NURBS: With Historical Perspective. San Francisco, CA: Morgan Kaufmann, 2001

[13] Lee S J, Park K S. Vibrations of Timoshenko beams with isogeometric approach. Applied Mathematical Modelling, 2013, 37(22): 9174–9190

[14] Hughes T J R, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39–41): 4135–4195

[15] de BoorC. A Practical Guide to Splines. New York, NY: Springer-Verlag, 1978

[16] Cottrell J A, Reali A, Bazilevs Y, Hughes T J R. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering, 2006, 195(41–43): 5257–5296

[17] Kheladi Z, Hamza-Cherif S M, Ghernaout M E A. Critical speeds analysis of spinning laminated composite shaft based on isogeometric analysis. International Journal for Computational Methods in Engineering Science and Mechanics, 2022, 23(6): 487–509

[18] Deng Q, Calo V M. A boundary penalization technique to remove outliers from isogeometric analysis on tensor-product meshes. Computer Methods in Applied Mechanics and Engineering, 2021, 383: 113907

[19] Hiemstra R R, Hughes T J R, Reali A, Schillinger D. Removal of spurious outlier frequencies and modes from isogeometric discretizations of second-and fourth-order problems in one, two, and three dimensions. Computer Methods in Applied Mechanics and Engineering, 2021, 387: 114115

[20] Manni C, Sande E, Speleers H. Application of optimal spline subspaces for the removal of spurious outliers in isogeometric discretizations. Computer Methods in Applied Mechanics and Engineering, 2022, 389: 114260

[21] Cheng L, Nicolas J. Free vibration analysis of a cylindrical shell-circular plate system with general coupling and various boundary conditions. Journal of Sound and Vibration, 1992, 155(2): 231–247

[22] Ertas B, John V. The influence of same-sign cross-coupled stiffness on rotordynamics. Journal of Vibration and Acoustics-Transactions of the ASME, 2007, 129(1): 24–31

[23] Cazzani A, Stochino F, Turco E. An analytical assessment of finite element and isogeometric analyses of the whole spectrum of Timoshenko beams. Journal of Applied Mathematics and Mechanics, 2016, 96(10): 1220–1244

[24] Genta G. Consistent matrices in rotor dynamic. Meccanica, 1985, 20(3): 235–248

[25] PieglLTillerW. The NURBS Book. Berlin: Springer Science & Business Media, 2012

[26] Yücel E, Saruhan H. Design optimization of rotor-bearing system considering critical speed using Taguchi method. Proceedings of the Institution of Mechanical Engineers Part E––Journal of Process Mechanical Engineering, 2017, 231(2): 138–146

Competing interests

The authors declare that they have no competing interests.

RIGHTS & PERMISSIONS

2024 Higher Education Press