Shrink-fitted joint behavior using three-dimensional solid finite elements in rotor dynamics with inclusion of stress-stiffening effect

Abstract

In this article, the contact behavior of a shrink-fitted joint in a rotor assembly at different spin speeds is studied, using the three-dimensional solid finite element approach. A custom frictionless contact model is proposed and extensively tested by means of simulation. The theory for solving pre-stressed damped eigenvalue problem is presented and implemented, and the results are presented in Campbell diagram form. Two sample rotor assemblies are studied. In addition, experimental modal analysis results of a shrink-fitted joint with various interferences are presented and numerically studied. The reduction of the interference of the shrink-fitted joint due to centrifugal forces at high-speed operation as well as the contact status are updated at each rotational speed step. The inclusion of stress-stiffening effect is studied in detail. The comparison of the Campbell diagram results with the results obtained using a commercial finite element software (Ansys) shows a good agreement.

Introduction

In high-speed machinery applications, rotor dynamic analysis is important part of the design progress. Often used one-dimensional (1D) beam element approach in rotor dynamics can represent only flexible shafts while the disks and impellers are modeled as rigid mass points. ¹ These 1D elements do not take into account the cross-section deformation. In certain applications, effects of deformable shaft cross section and flexible disks, and turbine blades are required. These features enable the inclusion of the stress-stiffening effect when the total stiffness of the rotor assembly is affected by the internal stress of the component. ² These aforementioned effects can be included when using two-dimensional (2D) axisymmetric harmonic or three-dimensional (3D) solid finite elements (FE).^3–5 The use of a 2D element is usually limited to axisymmetric structures and axisymmetric loads, whereas 3D elements can be used to describe, in practice, any rotating structure.

Use of 3D solid FEs in rotor dynamic analysis allows the inclusion of nonlinear contact behavior in the study of a rotor assembly. Modern rotor assemblies used in high-speed rotating machinery can include a large number of parts with various types of joints. The contact modeling and inclusion of frictional joint-based excitations are topics that have been challenging to include in the analysis when using the conventional beam element approach. Therefore, the use of a solid 3D element is the key to study contact-based phenomena. Additional phenomena such as centrifugal stiffening or stress stiffening provides further accuracy for the FE problem. Centrifugal stiffening is derived from the centrifugal load only, whereas stress stiffening is derived from all external loads taken into consideration.

General FE contact problems are discussed in Zhang and Wang, ⁶ where translational joint contact with friction force is studied in 2D, and in Nejati et al., ⁷ where internal friction of contact on a cracked region is studied. Many publications that consider the use of traditional beam or similar 1D elements such as Santos et al. ⁸ and Chen et al. ⁹ typically present a case-dependent solution for shrink-fitted sleeve–shaft contact. The solution is mainly about tuning the sleeve material properties, as based on the theory developed.

In the literature, Jafri ¹⁰ studied the sub-synchronous rotor dynamic instability caused by the shrink fit interface. His experimental studies showed that above the first critical speed (CS) of the rotor, strong unstable sub-synchronous vibration occurred, due to the shrink fit interface. Jafri ¹⁰ stated that the unstable sub-synchronous vibration originates from friction forces which are developed by the slippage in the shrink fit interface. These slippage-induced friction forces act as destabilizing cross-coupled moments while the rotor is operating above the first CS.

Chen et al. ¹¹ performed pre-stressed modal study of a hollow shaft assembly with shrink fits using Ansys. They used the native TARGE170 and CONTA174 contact elements of the software, and in order to verify the numerical results, experimental modal analysis (EMA) was performed. A similar study was presented in Chen, ¹² where an optimal equivalent direct model was proposed for manipulating the local stiffness in the contact region, by optimizing the elastic modulus of the material. Sasek ¹³ studied the sensitivity of eigenvalues of flexible rotors. He proposed a method of using individual subsystems for flexible disk and shaft and connecting these two subsystems using a global coupling force vector.

Other interesting types of contacts that are used in rotating machines are studied in Richardson et al. ¹⁴ and Qin et al. ¹⁵ In Richardson et al., ¹⁴ the bolt jointed Curvic coupling is under investigation and being studied using 3D FE method. The Curvic coupling is used especially in aero-engines for transmitting torque between rotor sections using a set of matching axial teeth having similar analogy than the hirth coupling. In Qin et al., ¹⁵ bolted disk–drum joints are studied using a formulated analytical model and the FE approach using Ansys software.

The effect of stress stiffening is shown as an important aspect to take into account in any FE-based frequency analysis. Studies of Kupnik et al. ¹⁶ show how significantly the natural frequencies of a gas turbine blade can be affected by inclusion of the stress-stiffening effect. Donadon et al. ¹⁷ present the results of change in harmonic response frequencies of an ultrasonic transducer while the stress-stiffening effect is included.

Still, there seems to be a complete lack of publications regarding the 3D solid element-based contact problem of a rotor assembly analysis in the rotor dynamics field. Based on previous studies of 3D solid FE rotor assemblies with shrink fits, the presence of commercial software seems to be dominant. Since the commercial software is designed to be able to solve the arbitrary contact problem, separate contact–target elements are used. The number of equations in rotor dynamic problems of rotor assemblies using finite solid elements is typically enormous. The use of separate contact detection element will only increase the computational effort of the solution. In this article, the shrink fit joint contact behavior is studied using 3D solid FEs in a MATLAB environment using a customized methodology. A direct node-to-node frictionless contact modeling approach is introduced in this article and used in the presented studies. The experimental modal analyses are performed to manufactured test shaft assemblies, and their dynamic behavior at high-speed operation is studied. Also, the contact region nonlinearity of a conical impeller assembly is studied. In addition, the effect of stress stiffening is studied in the aforementioned cases.

Method for modeling a pre-stressed rotor assembly

The methods for obtaining damped eigenvalues in an FE system with external forces using stress stiffening approach is presented. Contact detection without the use of separate contact–target elements and the contact force calculation procedure are proposed.

Constructing structural matrices

In this subsection, a short focused description of general construction of global structural mass and stiffness matrices are given. Global matrices are formed by assembling all constructed element matrices. Element mass and elastic stiffness matrices are given as follows^18,19

$${\mathbf{m}}_e=\underset{V}{\int }\rho {\mathbf{N}}^T\mathbf{N}\mathrm{dV}$$

(1)

$${\mathbf{k}}_e=\underset{V}{\int }{\mathbf{B}}^T\mathbf{DB}\mathrm{dV}$$

(2)

where $\mathbf{N}$ is the shape function matrix, $\rho$ is the material density, $\mathbf{B}$ is the strain matrix, and $\mathbf{D}$ is the matrix of material constants. The shape function matrix is given as follows

$$\mathbf{N}=\left[\begin{array}{cccc}{\mathbf{N}}_1\hfill & {\mathbf{N}}_2\hfill & \cdots \hfill & {\mathbf{N}}_n\hfill \end{array}\right]$$

(3)

where $n$ is the number of element shape functions. A submatrix ${\mathbf{N}}_i$ is given as follows

$${\mathbf{N}}_i=\left[\begin{array}{ccc}{N}_i\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_i\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & N_i\hfill \end{array}\right]$$

(4)

where $N_i$ is the shape function of node $i$ . The strain matrix is formed using the partial derivatives of the shape functions as follows

$$\mathbf{B}=\left[\begin{array}{ccc}\partial /\partial x\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \partial /\partial y\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \partial /\partial z\hfill \\ \multicolumn{1}{c}{\partial /\partial y}& \partial /\partial x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \partial /\partial z\hfill & \partial /\partial y\hfill \\ \multicolumn{1}{c}{\partial /\partial z}& 0\hfill & \partial /\partial x\hfill \end{array}\right]\mathbf{N}$$

(5)

The matrix of material constants for linear strain is given as follows

$$\mathbf{D}=\left[\begin{array}{cccccc}1-D_1\hfill & D_1\hfill & D_1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{D_1}& 1-D_1\hfill & D_1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{D_1}& D_1\hfill & 1-D_1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & D_2\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & D_2\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & D_2\hfill \end{array}\right]$$

(6)

where

$$D_1=-\frac{E\nu }{({\nu }^2-1)(1+\nu )}$$

(7)

$$D_2=\frac{E}{2(1+\nu )}$$

(8)

where $E$ is the material elastic modulus and $\nu$ is Poisson’s ratio. Additional steps for numerical integration and the use of isoparametric shape functions are given in Liu and Quek ¹⁸ and Logan. ¹⁹

Pre-stressed damped eigenvalue problem

A pre-stressed damped eigenvalue problem is used to study the rotor assembly with shrink fits. The equation of motion is as follows

$$\mathbf{M}\stackrel{··}{\mathbf{x}}+({\mathbf{C}}_e+{\mathbf{C}}_b+\mathrm{\Omega }\mathbf{G})\stackrel{·}{\mathbf{x}}+({\mathbf{K}}_e+{\mathbf{K}}_b+{\mathbf{K}}_c+{\mathbf{K}}_{\sigma })\mathbf{x}=0$$

(9)

In order to solve a damped eigenvalue problem, the following state space matrix A can be formed from equation (9)

$$\mathbf{A}=\left[\begin{array}{cc}0\hfill & \mathbf{I}\hfill \\ \multicolumn{1}{c}{-{\mathbf{M}}^{-1}{\mathbf{K}}_{\mathrm{tot}}}& -{\mathbf{M}}^{-1}({\mathbf{C}}_e+{\mathbf{C}}_b+\mathrm{\Omega }\mathbf{G})\hfill \end{array}\right]$$

(10)

where

$${\mathbf{K}}_{\mathrm{tot}}=({\mathbf{K}}_e+{\mathbf{K}}_b+{\mathbf{K}}_c+{\mathbf{K}}_{\sigma })$$

(11)

where $\mathbf{M}$ is the mass matrix, ${\mathbf{K}}_e$ and ${\mathbf{C}}_e$ are the unconstrained elastic stiffness and damping matrices, ${\mathbf{K}}_b$ and ${\mathbf{C}}_b$ are the bearing stiffness and damping matrices, correspondingly, ${\mathbf{K}}_c$ is the contact stiffness matrix, ${\mathbf{K}}_{\sigma }$ is the stress stiffness matrix, $\mathrm{\Omega }$ is the rotational velocity of the rotor, and $\mathbf{G}$ is the gyroscopic matrix. The solid FE theory that is used in this study is based on Liu and Quek ¹⁸ , Logan ¹⁹ and Zhou and Li. ²⁰ Due to circular shapes and usefulness of volumetric biasing of element size, the 10-node quadratic tetrahedron element type is used. The order of used shape functions is given in ANSYS, Inc. ²¹ Eigenvalues and eigenvectors are given by solving the eigenvalue problem ²²

$$\mathbf{Az}=\lambda \mathbf{z}$$

(12)

where $\mathbf{z}$ is the eigenvector and $\lambda$ is the eigenvalue of the particular eigenmode. These eigenvalues appear in a form of complex conjugate pairs, the ith eigenvalue pair is given

$${\lambda }_i={\alpha }_i\pm i{\beta }_i$$

(13)

where ${\alpha }_i$ and ${\beta }_i$ are the real and imaginary parts of eigenvalue $i$ , correspondingly and ${\beta }_i$ is the damped natural frequency of ith eigenmode.

In order to include the effect of closed contact in the eigenvalue problem, the contact region source–target nodes are constrained using node-to-node stiff spring connections. In eigenvalue analysis, matrix ${\mathbf{K}}_c$ is used to describe the surface-to-surface connections.

The stress stiffening matrix ${\mathbf{K}}_{\sigma }$ for a single element is as follows ²³

$${\mathbf{k}}_{\sigma }=\underset{V}{\int }{\mathbf{B}}_{ϵ}^T{\mathbf{S}}_0{\mathbf{B}}_{ϵ}\mathrm{dV}$$

(14)

where the matrix {\bf B}_{\epsilon} contains the spatial derivatives of the element shape functions

$$\mathbf{B}=\left[\begin{array}{cccc}{\mathbf{I}}_3\frac{\partial N_1}{\partial x}\hfill & {\mathbf{I}}_3\frac{\partial N_2}{\partial x}\hfill & \cdots \hfill & {\mathbf{I}}_3\frac{\partial N_n}{\partial x}\hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_3\frac{\partial N_1}{\partial y}}& {\mathbf{I}}_3\frac{\partial N_2}{\partial y}\hfill & \cdots \hfill & {\mathbf{I}}_3\frac{\partial N_n}{\partial y}\hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_3\frac{\partial N_1}{\partial z}}& {\mathbf{I}}_3\frac{\partial N_2}{\partial z}\hfill & \cdots \hfill & {\mathbf{I}}_3\frac{\partial N_n}{\partial z}\hfill \end{array}\right]$$

(15)

and the initial stress matrix ${\mathbf{S}}_0$ is

$${\mathbf{S}}_0=\left[\begin{array}{ccc}{\mathbf{I}}_3{\sigma }_x\hfill & {\mathbf{I}}_3{\tau }_{\mathrm{xy}}\hfill & {\mathbf{I}}_3{\tau }_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_3{\tau }_{\mathrm{yx}}}& {\mathbf{I}}_3{\sigma }_y\hfill & {\mathbf{I}}_3{\tau }_{\mathrm{yz}}\hfill \\ \multicolumn{1}{c}{{\mathbf{I}}_3{\tau }_{\mathrm{zx}}}& {\mathbf{I}}_3{\tau }_{\mathrm{zy}}\hfill & {\mathbf{I}}_3{\sigma }_z\hfill \end{array}\right]$$

(16)

where ${\mathbf{I}}_3$ is a 3 × 3 identity matrix.

The deformed element stress is calculated as follows ¹⁹

$$\{{\sigma }_x{\sigma }_y{\sigma }_z{\tau }_{\mathrm{xy}}{\tau }_{\mathrm{yz}}{\tau }_{\mathrm{zx}}{\}}^{\mathrm{T}}=\mathbf{DBu}$$

(17)

where $\{{\sigma }_x{\sigma }_y{\sigma }_z\}$ are the element normal stress components, $\{{\tau }_{\mathrm{xy}}{\tau }_{\mathrm{yz}}{\tau }_{\mathrm{zx}}\}$ are the element shear stress components, $\mathbf{D}$ is the elasticity matrix, and $\mathbf{B}$ is the element strain matrix. The element nodal displacement vector $\mathbf{u}$ is calculated as follows

$$\mathbf{u}={\mathbf{K}}_{e,\mathrm{CYL}}^{-1}({\mathbf{\Omega }}^2{\mathbf{F}}_{\mathbf{\Omega }}+{\mathbf{F}}_c)$$

(18)

where ${\mathbf{K}}_{e,\mathrm{CYL}}$ is the elastic stiffness matrix with cylindrical constraints allowing radial and axial deformations but restricting tangential deformations in accordance with the axis of rotation, and ${\mathbf{F}}_c$ is the contact force vector. The centrifugal force vector ${\mathbf{F}}_{\mathrm{\Omega }}$ is calculated as follows ²³

$${\mathbf{F}}_{\mathrm{\Omega }}=\rho \underset{V}{\int }{\mathbf{N}}^T{\mathbf{\Omega }}^T\mathbf{\Omega }x\mathrm{dV}$$

(19)

where $\rho$ is the material density, $\mathbf{N}$ is the element shape function matrix, and $\mathbf{x}$ is the vector of nodal coordinates. The matrix $\mathbf{\Omega }$ indicating the axis of rotation is as follows

$$\mathbf{\Omega }\mathbf{=}\mathbf{\left[}\begin{array}{ccc}0\hfill & \mathbf{-}{\mathbf{\omega }}_{\mathbf{z}}\hfill & {\mathbf{\omega }}_{\mathbf{y}}\hfill \\ \multicolumn{1}{c}{{\mathbf{\omega }}_{\mathbf{z}}}& 0\hfill & \mathbf{-}{\mathbf{\omega }}_{\mathbf{x}}\hfill \\ \multicolumn{1}{c}{\mathbf{-}{\mathbf{\omega }}_{\mathbf{y}}}& {\mathbf{\omega }}_{\mathbf{x}}\hfill & 0\hfill \end{array}\mathbf{\right]}$$

(20)

where ${\omega }_x{\omega }_y{\omega }_z$ are the normalized angular velocity components, according to the global coordinate system, $|\{{\omega }_x{\omega }_y{\omega }_z\}|=1$ .

Contact modeling

The interference caused by shrink fit in a rotor assembly can be modeled as either geometrical or numerical interference between shaft and sleeve bodies. The numerical approach is more convenient when, for example, the optimizing routine is used to find a suitable interference fit. When using numerical interference, the bodies are modeled without geometrical interference, and there is no penetration or gap in the contact region. In addition to other possible external forces, the contact force is applied between surfaces of contact bodies in order to yield the desired radial gap to the contact region. The radial gap is then equal to or greater than the designed radial interference, thereby allowing contact to be closed or open.

In this article, elementless frictionless contact modeling is proposed, since it enables the use of less FEs resulting in reasonably lower computational effort. The main core of elementless contact detection is the fact that the FE mesh on the contact region of both the source and target bodies, as illustrated in Figure 1, is refined so that direct node-to-node penetration detecting can be utilized. Thus, there is no need for separate contact elements. Friction forces on the contact region are neglectable due to the lack of any external axial or torsional forces allowing usage of the frictionless contact model. Studied FE problems are meshed with converging node distribution in the contact region. This means that the distance between source–target nodes is zero before applying the radial interference, if no external forces are present. This proposed node-to-node contact between the surfaces of two bodies is flexible-to-flexible contact by default.

Figure 1.

Source and target node i of sleeve and shaft body.

The contact model introduced in this article uses a fixed number of evaluations. At this point, it is assumed that no sliding or slipping occurs in the contact region. The contact detection is done using the cylindrical coordinate system. The deformations are mapped from the Cartesian coordinate system into the cylindrical coordinate system. After the contact force vector is calculated, it is mapped back into the Cartesian coordinate system. The first evaluation according to equation (21) takes into account all the external forces that, in these studied cases, represent centrifugal force only. The deformation vector for the first evaluation is calculated as follows

$${\mathbf{u}}_0={\mathbf{K}}_{e,\mathrm{CYL}}^{-1}\left({\mathrm{\Omega }}^2{\mathbf{F}}_{\mathrm{\Omega }}\right)$$

(21)

The radial distance between source–target nodes $\mathbf{\Delta }{\mathbf{u}\prime }_0$ is given as follows

$$\mathbf{\Delta }{\mathbf{u}\prime }_0={\mathbf{u}\prime }_{0,s}-{\mathbf{u}\prime }_{0,t}$$

(22)

where the superscript ′ describes the cylindrical coordinate system and the subscripts $s$ and $t$ correspond to the source and target as seen in Figure 1. The vector $\mathbf{\Delta }{\mathbf{u}}_0$ is converted from the global Cartesian coordinate system into the cylindrical coordinate system. The vector of logical multipliers $\gamma {\prime }_0$ denotes which contact node pairs are opened and which are closed and can be defined as follows

$$\gamma {\prime }_{0,u_i}=\{\begin{array}{cc}1,\hfill & \mathrm{if}\mathrm{\Delta }{u}_{0,i}\prime <{\delta }_{\mathrm{rad}}\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{if}\mathrm{\Delta }{u}_{0,i}\prime ⩾{\delta }_{\mathrm{rad}}\hfill \end{array}$$

(23)

where $\gamma {\prime }_{0,u_i}$ is the multiplier of radial component of source–target node pair $i$ and ${\delta }_{\mathrm{rad}}$ is the radial interference of the shrink fit joint. The second evaluation step uses the inclusion of centrifugal force and trial contact force. Two known variables, centrifugal force and trial contact force, help in determining the unknown contact stiffness. By solving the contact stiffness, the unknown final contact force can be calculated. The deformation vector for the second evaluation is calculated as follows

$${\mathbf{u}}_{\mathrm{trial}}={\mathbf{K}}_{e,\mathrm{CYL}}^{-1}({\mathrm{\Omega }}^2{\mathbf{F}}_{\mathrm{\Omega }}+{\mathbf{F}}_{c,\mathrm{trial}})$$

(24)

where ${\mathbf{F}\prime }_{c,trial}$ is the vector of trial contact force in the cylindrical coordinate system given as

$${\mathbf{F}\prime }_{c,trial}=k_0\mathbf{\Delta }{\mathbf{u}}_0·\gamma {\prime }_0$$

(25)

where $k_0$ is the initial value for the contact stiffness. The final contact force vector is calculated as follows

$${\mathbf{F}\prime }_c={\mathbf{k}}_{\mathbf{f}}{k}_0·\mathbf{\Delta }{{\mathbf{u}}^{\prime }}_0·{{\gamma }^{\prime }}_0$$

(26)

where ${\mathbf{k}}_{\mathbf{f}}$ is the vector of multipliers for initial contact stiffness. When using the linear material model, the multiplier for source–target node pair $i$ is given by

$${\delta }_{\mathrm{rad}}-\mathrm{\Delta }u{\prime }_{0,i}-k_{f,i}(\mathrm{\Delta }{u\prime }_{\mathrm{trial},i}-\mathrm{\Delta }{u\prime }_{0,i})\gamma {\prime }_{0,u_i}=0$$

(27)

If $\gamma {\prime }_{0,u_i}=1$ , then

$$k_{f,i}=\frac{{\delta }_{\mathrm{rad}}-\mathrm{\Delta }{u\prime }_{0,i}}{\mathrm{\Delta }{u\prime }_{\mathrm{trial},i}-\mathrm{\Delta }{u\prime }_{0,i}}$$

(28)

The vector $\mathbf{\Delta }\mathbf{u}{\prime }_{\mathrm{trial}}$ calculated as based on the trial evaluation is given as follows

$$\mathbf{\Delta }{\mathbf{u}\prime }_{\mathrm{trial}}=\mathbf{u}{\prime }_{\mathrm{trial},s}-\mathbf{u}{\prime }_{\mathrm{trial},t}$$

(29)

Experimental results from the test shaft assembly

The test shaft assembly under study consists of a steel shaft and shrink-fitted aluminum sleeve, as illustrated in Figure 2. Certain material properties of individual components are verified by measuring mass properties and the lowest free–free natural bending mode (BM) frequencies (shafts only). For shrink-fitted joint testing purposes, three sets of test assemblies are manufactured with different radial interferences and measured by means of EMA. A scanning laser doppler vibrometer (SLDV) is used, which provides contactless measurement and does not require sensors that would add mass on the test pieces. The experimental results are shown in Table 1. Based on the measured first free–free natural BM frequencies, the radial interference in the shrink fit joint does not exert a significant effect on the natural frequencies of the test assemblies.

Figure 2.

Test shaft assembly during experimental modal analysis.

Table 1.

Experimental results of different radial interferences of test shaft assemblies.

Test shaft assembly	Radial interference, δ rad (µm)	1. BM (Hz)	2. BM (Hz)
1	3.0	1105.9	3084.0
2	9.0	1108.2	3098.4
3	14	1108.6	3095.7
Average		1107.6	3092.7
Difference (max/min) a		+0.09%/−0.15%	+0.18%/−0.28%

BM: bending mode.

a

Max/min value compared against the corresponding averaged value.

Component-level verification of certain material properties was performed before assembling the shrink fits. Table 2 shows the material parameters verified and applied in the numerical studies. The density of both materials used is verified by measuring the masses of all test pieces and by comparing the averaged results against known volume. The known volume is taken from the computer-aided design (CAD) model. Manufactured bodies will have marginally different dimensions due to the tolerances; however, the change in volume between test pieces was considered insignificant. The same conclusion was made while comparing the volumes of the CAD body and FE mesh. The elastic modulus of S235 is verified by measuring the first natural bending frequencies of all shaft parts. The averaged results as seen in Table 1 are set as a reference. The elastic modulus of S235 is found using the iteration loop in FE software (Ansys). The elastic modulus of 6061-T6 aluminum was not verified. The value used for elastic modulus of 6061-T6 is taken from the material data sheet. Also, Poisson’s ratios are based on material data sheets for both materials.

Table 2.

Material parameters used with numerical studies of test shaft assembly.

Part/material	Density (kg/m3)	Elastic modulus (GPa)	Poisson’s ratio (−)
Shaft/S235	7812 a	210 a	0.3
Sleeve/Al 6061-T6	2745 a	68.9	0.33

a

Verified parameter.

Numerical results of test shaft assembly

Based on the measured results of manufactured test shaft assemblies, the shrink fit interference seems to have a minor role in the free–free natural bending frequencies. Therefore, the case of 3 µm interference was selected for numerical studies. The studied test shaft assembly spin speed range is from 0 to 30,000 r/min with 1000 r/min increment. Another reason for selecting this particular test case is the fact that the shrink fit joint becomes loose as the rotation speed reaches 23,000 r/min, and as a result the maximum cannot be reached.

The numerical results are obtained using the proposed method, and results from commercial software Ansys Workbench 18 are obtained for comparison purposes. The contact treatment used in the MATLAB environment using the proposed method is based on the theory introduced. Results are shown as Campbell diagrams and tabulated as natural frequencies. The effect of stress stiffening is explored by studying the system with and without its influence. Thermal effects are not included in any of the presented studies.

A comparison of free–free natural BMs is shown in Table 3. The case of 3 µm interference is studied. In the case of the method proposed, results using frictionless contact with and without the initial stress stiffening, due to shrink fit interference, are given. Results obtained using both the proposed method and Ansys solution with frictional contact are compared against the averaged results of the measured test shaft assemblies. The averaging of the measured results is justified, since the effect of radial interference in the shrink fit joint does not significantly affect the natural free–free bending frequencies of the measured test shaft assemblies. The proposed method yields fairly accurate results, though Ansys reference results yield greater error compared to the measured results. Since the material properties used are the same with both modeling approaches, the contact modeling is logically the main reason for differing results. The material properties used in the numerical studies are given in Table 2.

Table 3.

Comparison of test shaft assembly free–free bending frequencies.

	1. BM (Hz)/error (%)	2. BM (Hz)/error (%)
Measured (averaged)	1107.6/–	3092.7/–
Proposed method without stress stiffening	1113.4/+0.52	3147.0/+1.76
Proposed method with stress stiffening	1113.3/+0.51	3147.0/+1.76
Ansys	1138.3/+2.77	3300.8/+6.73

BM: bending mode.

In addition, the test shaft assembly cases 2 and 3 according to Table 1 are studied. Free–free bending frequency results are calculated with proposed method and Ansys. Between the different radial interferences, the numerical results are the same. When the stress-stiffening effect was included with the proposed method, the bending frequencies increased only by 0.02%, due to increased joint interference.

The FE-mesh of the test shaft assembly used with the proposed method is presented in Figure 3. The Ansys Workbench is used for generating FE-mesh. The mesh consist of 36,838 ten-node quadratic tetrahedron elements. The mesh size at the contact zone is 2.5 mm. The shaft diameter at the contact zone is 20 mm. By using a mesh imported from Ansys, possible differences of model discretization are avoided, and therefore differences between two modeling approaches are minimized. The length of the shaft is 260 mm and diameter of the sleeve is 80 mm. Two simple identical bearings, essentially spring-damper elements, are located at both ends of the shaft and connected to the surfaces of end planes of the shaft. The used bearing parameters are presented in Table 4.

Figure 3.

FE-mesh of test shaft assembly.

Table 4.

Bearing data used in the test shaft assembly.

Bearing	Stiffness		Damping
	kyy (N/m)	kzz (N/m)	cyy (N s/m)	czz (N s/m)
1	1 × 108	1 × 108	1 × 103	1 × 103
2	1 × 108	1 × 108	1 × 103	1 × 103

The Campbell diagram generated using the proposed method without the stress-stiffening effect is presented in Figure 4. It is known that the minor radial interference of 3 µm in the shrink fit joint will not hold the sleeve–shaft contact closed during the full speed range. The threshold is 23,000 r/min as seen in Figure 4. After 23,000 r/min, the contact is fully open and the resulting eigenvalues would be those of two distinct parts (shaft and sleeve) that are not comparable to the eigenvalues of the assembly. Therefore, above the threshold rotational speed, the natural frequencies are not gathered in the Campbell diagram. Also, when the contact opens, six additional rigid body modes can be obtained. The backward whirl (BW) and forward whirl (FW) mode frequencies at the threshold speed of 23,000 r/min are marked in Figures 4–7.

Figure 4.

Campbell diagram without stress stiffening of the test shaft assembly, zero speed bending mode (BM) frequency, backward whirl (BW) and forward whirl (FW) frequencies marked.

Figure 5.

Proposed method Campbell diagram with stress stiffening of the test shaft assembly.

Figure 6.

Proposed method Campbell diagram of the test shaft assembly using fixed contact, critical speeds (CS), and corresponding whirling modes marked.

Figure 7.

Ansys Campbell diagram of the test shaft assembly using frictional contact.

In Figure 5, the stress-stiffening effect is included. The stiffening effect slightly increases the frequencies of the first forward and BW modes. The Campbell diagram calculated with fixed contact, and therefore without initial contact interference nor the stress-stiffening effect is presented in Figure 6 for comparison purposes. The contact will not open at high speed due to fixed contact between shaft and sleeve. The CSs and corresponding modes are marked in Figures 6 and 7.

The Ansys Campbell diagram of the test shaft assembly is presented in Figure 7. The frictional contact with initial radial interference in the shrink fit joint is applied. As it can be clearly seen in Figure 7, Ansys does not take into account the contact opening. The results between the proposed method with fixed contact and Ansys results with frictional contact behave very similarly.

Zero speed-supported test shaft assembly frequencies of the first BM are the same using frictionless and fixed contact with the proposed method. Ansys predicts this first BM to be 9 Hz greater at zero speed. The 23,000 r/min speed was considered as one comparison point, since beyond this speed the contact treatment detects the contact fully opened. For this reason, the rest of the Campbell diagram cannot be solved when using frictionless contact with the proposed method. At 23,000 r/min, the proposed method results are consistent: the stress-stiffening effect increases whirling frequencies by 7 Hz, whereas Ansys predicts over 9 Hz greater whirling frequencies. While comparing the frictionless method results against the Ansys results, the latter ones seems constantly to be 9 Hz greater between 0 and 23,000 r/min. The results available at the upper limit of the studied speed range (30,000 r/min) are compiled in Table 5 for comparison purposes. At 30,000 r/min, only the results available using the proposed method are the results using the fixed contact model. It can be observed from the table that whirl frequencies obtained from the Ansys model are still over 8 Hz greater than the ones predicted using the proposed method. The reason for this can be in the contact modeling that Ansys uses. Similar behavior—the greater free–free bending frequencies—is shown in Table 3. Critical BW and FW speeds of the Ansys model with frictional contact are a bit over 500 r/min greater than the CSs predicted using the proposed method with fixed contact. Frictionless contact model with the proposed method does not reach the first CSs.

Table 5.

Whirling mode frequencies at 30,000 r/min of the test shaft assembly.

	At 30,000 r/min
	1. BW (Hz)	1. FW (Hz)
Proposed method with fixed contact	431.1	476.6
Ansys	440.9	484.8

BW: backward whirl; FW: forward whirl.

Inclusion of the stress-stiffening effect will increase the mode frequencies by 7 Hz at 23,000 r/min. Also, the frequency difference due to the gyroscopic effect between BW and FW modes at 23,000 r/min using the proposed method is practically the same: average of 34.65 Hz ± 0.44%. The contact modeling method and stress stiffening clearly increase the overall stiffness of the structure.

Conical impeller assembly

The construction of the conical impeller shrink-fitted on the shaft is studied. The FE-mesh of the conical impeller assembly is presented in Figure 8. The mesh consists of 24,766 ten-node quadratic tetrahedron elements. The mesh size at the contact zone is 5.0 mm. The shaft diameter at the contact zone is 60 mm. Length of the shaft is 800 mm and the maximum diameter of the impeller wheel is 300 mm. The conical impeller of the radial kinetic compressor wheel is modeled without the blades. Only the impeller hub is included in the FE-mesh. By doing so, the size of the FE-problem can be reduced and the lower number of the lowest eigenmodes can be utilized by neglecting the low-frequency local BMs of the impeller blades. The material properties used in this study are presented in Table 6.

Figure 8.

FE-mesh of conical impeller assembly.

Table 6.

Material parameters used in numerical studies of conical impeller assembly.

Part/material	Density (kg/m3)	Elastic modulus (GPa)	Poisson’s ratio (−)
Shaft/steel	7850	200	0.3
Impeller/aluminum	2810	71.7	0.33

Two simple identical bearings, essentially spring-damper elements, are located at both ends of the shaft and are connected to the surfaces of end planes of the shaft. The stiffness and damping parameters used are the same as in the test shaft assembly example presented in Table 4. The studied speed range is from 0 to 30,000 r/min with 1000 r/min increment. The radial interference in the impeller–shaft shrink fit is 50 µm.

This particular simplified conical impeller–shaft assembly is selected because the impeller–shaft shrink-fitted joint will partially open as a function of spin speed. Partially opened contact region is illustrated in Figure 9 while the deformation is amplified for visualization purposes. The contact area will decrease while operating speed increases. This will cause a reduction of bending stiffness in this assembly.

Figure 9.

Visualization of contact region partial opening between the conical impeller and shaft.

Campbell diagram of the conical impeller assembly without stress-stiffening effect is presented in Figure 10. The effect of partially opened contact is present in the Campbell diagram as the overall stiffness of the assembly is starting to reduce at 15,000 r/min. A Campbell diagram of the conical impeller assembly with stress-stiffening effect is presented in Figure 11. The stress-stiffening effect will not yield a significant stiffening effect on the BMs.

Figure 10.

Proposed method: the Campbell diagram without stress stiffening of the conical impeller assembly.

Figure 11.

Proposed method: the Campbell diagram with stress stiffening of the conical impeller assembly.

The Campbell diagram of the conical impeller assembly with fixed contact is presented in Figure 12 for comparison purposes. The initial interference in the shrink fit joint and the stress-stiffening effect are neglected. Since the contact is fixed, the effect of stiffness reduction due to partially opened contact is not present in the Campbell diagram. The Campbell diagram generated in Ansys is presented in Figure 13. Both the fixed joint results in Figure 11 and the Ansys results behave very similarly without having any indication of reduced bending stiffness due to partially opened contact.

Figure 12.

Proposed method: the Campbell diagram of the conical impeller assembly using fixed contact.

Figure 13.

Ansys Campbell diagram of the conical impeller assembly using frictional contact.

With the conical impeller assembly, the first two BMs are visible in the Campbell diagrams (Figures 10–13). The first BW and FW CS of Ansys are consistently 50 r/min greater than the ones with the proposed method. Even so, the results obtained using the proposed method with the stress-stiffening effect are actually the same as the Ansys model frequencies. Due to contact treatment applied on every speed step while using the proposed method, the contact opening is taken into account. Due to reduction in the contact surface area at high speeds, the dip in the mode lines in the Campbell diagrams (Figures 10 and 11) is present. Only the models with frictionless contact using the proposed method experience the second BW mode in the speed range studied. The stress-stiffening effect increases the second BW CS by 360 r/min.

The reason for a rather large reduction in whirl mode frequencies of the second BM using frictionless contact is explicable due to the BM shapes of the supported structure. The two lowest zero-speed BMs of the supported conical impeller assembly are illustrated in Figure 14 using Paraview visualization software. The node points of the first BM are at bearing locations. In the case of the second BM, the third node appears roughly at the center of the structure. Because the impeller–shaft joint is located roughly one quarter of the assembly length away from the other end, the joint is roughly located at the position of highest amplitude of the second BM. Also the bending deformation of the shaft is at the highest in the location of the impeller and the decreasing bending stiffness has a significant effect in this particular mode.

Figure 14.

Conical impeller assembly lowest supported BMs at zero speed: (a) first supported BM at zero speed and (b) second supported BM at zero speed.

The zero-speed-supported conical impeller assembly frequencies are all very consistent. In Table 7, the whirling mode frequencies at 30,000 r/min are compiled. At 30,000 r/min, the first whirling mode frequencies are roughly the same, while the proposed method with stress stiffening yields greatest frequencies. The second whirling modes solved using frictionless contact with the proposed method are no longer comparable to the Ansys results, due to partially opened contact and, therefore, the reduced overall bending stiffness of the conical impeller assembly. An interesting observation is that the Ansys model with frictional contact behaves very similar to the model using the proposed method with fixed contact.

Table 7.

Whirling mode frequencies at 30,000 r/min of the conical impeller assembly.

	At 30,000 r/min
	1. BW (Hz)	1. FW (Hz)
Proposed method without stress stiffening	161.1	171.1
Proposed method with stress stiffening	168.6	176.8
Proposed method with fixed contact	165.5	173.3
Ansys	166.4	173.9
	At 30,000 r/min
	2. BW (Hz)	2. FW (Hz)
Proposed method without stress stiffening	484.7	511.7
Proposed method with stress stiffening	491.4	519.9
Proposed method with fixed contact	505.0	530.0
Ansys	508.0	533.4

BW: backward whirl; FW: forward whirl.

The partially opened contact between conical impeller–shaft will reduce the whirling mode frequencies. When comparing the second BM frequencies, the backward mode will reduce by average of 18.5 Hz and the forward mode by average of 15.9 Hz. The stress-stiffening effect will increase the stiffness of conical impeller and shaft in proximity of the contact region. Even so, the total stiffness of this structure is dominated by the robustness of the conical impeller–shaft joint. Since the contact is partially opened at full speed, the total stiffness of the structure is reduced regardless of the stress-stiffening effect. By comparing the forward and backward mode frequency difference at 30,000 r/min, the gyroscopic effect will also increase due to contact opening by an average of 10.1%.

Accuracy of contact modeling

In this section, the constraint error at the contact region achieved using the proposed contact modeling method is presented. The magnitude of the constraint error is relevant in order to achieve converged results. For comparison purposes, the results of the constraint error at the contact region calculated with Ansys are presented. Figure 15 introduces the cylindrical constraints that are used as boundary conditions for the contact treatment with the proposed method. The displacements of the assembly according to Figure 15 are constrained as follows: radial deformation is allowed while tangential deformation is restricted. Axial rigid body movement is restricted, but axial deformations are allowed. In Ansys, a cylindrical constraint that works the same way is applied.

Figure 15.

Cylindrical constraints applied: red dot TX = 0, green dots TY = 0, blue dots TZ = 0.

For a reference solution, the test shaft assembly is selected. Rotational speed is set to be zero to ensure that no external loads are present, and thus the contact region is fully closed. The constraint error with the proposed method at the contact region is presented in Figure 16.

Figure 16.

Constraint error at the contact region using proposed method.

The proposed contact method is based on the use of contact stiffness coefficient. Therefore, it resembles the Penalty function constraint equation method. ²⁴ With this method, it is typical that the solution is not exact because of the penalty weight, essentially a spring coefficient, is used. The magnitude of the cells in the global elastic stiffness matrix and the penalty weight used can be used to determine the amplitude of the maximum error in the constraint. ²⁴ As seen in Figure 16, the solution is not very focused, although the mean error is relatively small. The solution required two iterations. The cylindrical constraint utilizes approach with stiff springs that are connected to the ground. Therefore, the cylindrical constraint may also contribute a minor constraint error.

For comparison purposes, the constraint error at the contact region is calculated using Ansys. The Augmented Lagrange contact formulation is used. The constraint error in this case is presented in Figure 17.

Figure 17.

Constraint error at the contact region using Augmented Lagrange method with Ansys.

As seen in Figure 17, the constraint error while using the Augmented Lagrange method is relatively focused, although it contains minor error. The solution required four iterations. The results of the constraint error presented in this section will give a good general understanding how different constraint method work, how many iterations they require, and how accurate they are.

Conclusion

In this article, a custom-contact model is developed and utilized for studying contact behavior at high-speed operation. The results are presented and compared to the results obtained using commercial software (Ansys). The inclusion of the stress-stiffening effect is also studied extensively. Two example cases are studied: an experimentally verified test shaft assembly with shrink fit interference joint and a conical impeller assembly. With the test shaft assembly, a number of different radial interferences are studied. Based on the results, it was noticed that the amount of radial interference in the shrink-fitted joint will not in fact change the free–free natural bending frequencies of the assembly noticeably because of the low pre-stress state. The high-speed operation of shrink-fitted sleeve was studied with the custom contact model to the point when the contact will be fully opened. At greater speeds, the solution of eigenvalue problem would yield eigenvalues of individual bodies instead of the assembly.

With the conical impeller assembly, the partial contact opening was studied and detected using a custom frictionless contact model. The same conditions are studied also using Ansys software for comparison purposes. The custom contact model proposed will update the contact status on every speed step. The results generated using this custom contact model differ moderately from the ones solved using Ansys. It seems the contact status is not updated in Ansys during the regular modal analysis where the rotational speed is used as a variable forming eventually the ingredients for generating Campbell diagram. Based on the results presented in this article, the stress-stiffening effect should be included into eigenvalue calculation particularly in cases, where relatively heavy bodies attached on the shaft using interference fittings are rotating high-speeds.

The constraint error results using the proposed method and Ansys frictional contact, as part of convergence analysis, are presented and discussed. Future studies in contact modeling area would include a transient study of frictional contact taking into account the contact opening and the friction forces as a possible source of instability.