A Numerical Simulation of Elastoplastic Contact Analysis of Compressor by Overspeeding

Abstract

This study analyzes the stress distribution of 3D elastoplastic contact problems by using the FE parametric quadratic programming (PQP) method derived from a 3D FE model based on parametric variational principle (PVP). We numerically analyze a 24-blade compressor by combining centrifugal loading with interference-fit one. To accelerate computation, calculation is simplified by structural modeling via multisubstructuring, aiming to deal with FE-simulated computer aided design (CAD) conveniently. We then analyze the relationships between the maximum residual stresses of the compressor posterior to prestressing and overspeed rpms, and we also study the distribution and magnitude of the contact stresses of the compressor in working conditions after overspeed prestressing. Moreover, we thoroughly discuss the distribution and magnitude of the contact stresses of shaft-shaft sleeve-impeller in working conditions. Relative displacement can be prevented and contact stress can be kept uniform due to the nonuniform initial amount of interference in overspeed prestressing. This paper summarizes the FEM simulation results and provides reference data for improving the design and processing of compressor impellers, indicating that overspeed is indispensable in manufacture.

1. Introduction

Turbocharger, which is supported by two inboard-mounted full-floating shaft sleeves, is a device in which a radial inflow turbine is turned by the engine exhaust gas-originating heat energy with a centrifugal compressor on the same shaft. As a result, the combustion efficiency of engine is boosted owing to the pressurized inflow air provided by compressor [1]. Accordingly, nearly 75% of the total output of an engine stems from turbocharging, thus rendering its running stability prerequisite for other components.

Centrifugal compressor, which is the most available high-speed turbomachinery [2], is featured in high rotational speed, centrifugal force, and stability. Exhaust flow from engine is directed over turbine blades to drive the turning of compressor and shaft. Accelerating the turning of compressor impeller boosts engine more effectively. Hence, impeller is mainly subjected to centrifugal forces as a highly stressed component due to the high-speed rotation in operation. The roots of impeller blades are prone to plastic deformations because impeller is commonly made from forged aluminum alloy.

It is crucial to prestress impeller by overspeed for performance optimization. Impeller material is strengthened by plastic deformation (Figure 1), that is, yield limit increase, to expand elastic range in case of forced deformation, to raise yield point and carrying capacity, and to realize self-strengthening. Besides, the size of the area undergoing forced deformation cannot be recovered when overspeed ceases owing to the restraint of peripheral elastic area that leads to residual stress and deformation. Residual compressive stress, which is brought about in the area undergoing plastic deformation in overspeed prestressing, lowers the stress at normal working speed. Meanwhile, overspeed inhibits the extension of microscopic defects and substantially extends fatigue life by passivating the defects during heat treatment or processing. The method is described below.

Figure 1:

Stress-strain curve metal under loading, unloading, and reloading.

Impeller is rotated at a higher speed than the rated one before assembling, which thus increases plastic deformations (point B in Figure 1). When the impeller is unloaded and stopped, the stress-strain characteristic curve returns to point C along the straight line parallel to its elastic section rather than the coordinates origin along the constitutive relation curve, giving rise to residual stresses and deformations in the inner portion. During the impeller assembling, the stress-strain curve rises from point C instead of the coordinates origin to point D with aggravating residual deformation corresponding to that at the rated speed along the straight line CB. Obviously, the overspeed-prestressed impeller generates residual stress and deformation adjoining to the inner radius, thus optimizing the compressive stress state in the inner portion. During service, the stress state at lower speeds is superimposed to give a net and more favourable one, so that the impeller may function at higher working speeds compared with that induces yielding onset at the inner radius.

After overspeed prestressing, the impeller is mounted onto the shaft via interference fit, which burdens a complicated frictional contact problem on the interface of mutually connected parts at high rotational speed. Elastoplastic contact problems are tricky because their material and geometric nonlinearities are actually nonlinear for the moving boundary and friction along the contact surface. Several alike problems have recently been reported [3–8]; the numerical analysis methods for which, such as parametric variational principle (PVP) and the corresponding parametric quadratic programming (PQP), have also been proved successful. Zhong et al. [9, 10] developed a PVP to analyze and reduce plane contact problems to parametric programming problems by the finite element method and to linear complementarity problems finally [11–17]. This method outweighs conventional ones because of the eliminated penalty factors and the facile solutions in the absence of tedious iterations (e.g., general incremental iteration).

Three-dimensional elastoplastic contact problems with friction have been numerically treated [13] by being considered as a special case that constrain-minimizing either complementary or total potential energy of elastoplastic structures, and by eliminating penalty factors via special protocols. The conditions avoiding nonphysical ray solution and defining a unique solution map are given. Furthermore, a combined PQP and iterative algorithm are introduced to simplify the determination of slip direction, thus linearizing Coulomb's friction law and increasing the solution accuracy. Notably, the method described herein can be utilized for multistep loading incremental problems. The derivation is shown in a one-step form for clarification, and the formulations corresponding to multistep problems can be easily derived.

A 24-blade locomotive-type turbocharger compressor combining interference fit and centrifugal loadings is utilized in the PQP method for numerical modelling. The interference fit is analyzed based on PVP, and complex impeller geometrical shapes are simulated precisely with accurate numerical computation by the multilevel and multibranch substructuring of FE. The influences of fit tolerance, rotational speed (centrifugal force), and the wall thickness of shaft sleeve on the distribution of contact stress as well as structure are also thoroughly discussed.

2. Computational Theory

2.1. Governing Equations

Considering the orthotropic behaviors, the governing equations for 3D contact analysis can be summarized as follows.

Equilibrium equation:

$d{\sigma }_{ij,j}+d{b}_i=\mathrm{0}.$

(1)

Continuity equations:

$d{\epsilon }_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left(d{u}_{i,j}+d{u}_{j,i}\right).$

(2)

Constitutive equations:

$\begin{array}{c}d{\sigma }_{ij}=D_{ijkl}\left(d{\epsilon }_{kl}-d{\epsilon }_{kl}^p\right),\hspace{1em}f\left({\sigma }_{ij},{\epsilon }_{ij}^p,\kappa \right)\le \mathrm{0},\\ \hspace{1em}d{\epsilon }_{kl}^p=\frac{\partial g}{\partial {\sigma }_{kl}}\lambda ,\hspace{1em}\lambda \{\begin{array}{cc}\ge \mathrm{0},\hfill & f=\mathrm{0},\hfill \\ =\mathrm{0},\hfill & f<\mathrm{0}.\hfill \end{array}\end{array}$

(3)

Boundary conditions:

$\begin{array}{c}{n}_{i}d{\sigma }_{ij}=d{\stackrel{-}{p}}_j,\hspace{1em}\mathrm{on}\mathrm{\hspace{0.17em}\hspace{0.17em}}{S}_p;\\ d{u}_i=d{\stackrel{-}{u}}_i,\hspace{1em}\mathrm{on}\mathrm{\hspace{0.17em}\hspace{0.17em}}{S}_u;\\ {\epsilon }_n=d{u}_n^{\left(\mathrm{1}\right)}-d{u}_n^{\left(\mathrm{2}\right)}+{\delta }^{*}\ge \mathrm{0},\hspace{1em}\mathrm{on}\mathrm{\hspace{0.17em}\hspace{0.17em}}{S}_c;\\ {\text{if}\mathrm{\hspace{0.17em}\hspace{0.17em}}\epsilon }_n>\mathrm{0}⟹p_n=p_{{\tau }_{\mathrm{1}}}=p_{{\tau }_{\mathrm{2}}}=\mathrm{0};\\ {\text{if}\mathrm{\hspace{0.17em}\hspace{0.17em}}\epsilon }_n=\mathrm{0}⟹p_n\le \mathrm{0},\hspace{1em}{\left[p_{{\tau }_{\mathrm{1}}}^{\mathrm{2}}+{\beta }^{\mathrm{2}}{p}_{{\tau }_{\mathrm{2}}}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}+{\mu }_{\mathrm{1}}{p}_n\le \mathrm{0};\end{array}$

(4)

if

$\begin{array}{c}{\left[p_{{\tau }_{\mathrm{1}}}^{\mathrm{2}}+{\beta }^{\mathrm{2}}{p}_{{\tau }_{\mathrm{2}}}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}+{\mu }_{\mathrm{1}}{p}_n<\mathrm{0}⟹\left|\mathrm{\Delta }{u}_{\tau }\right|=\mathrm{0},\\ {\left[p_{{\tau }_{\mathrm{1}}}^{\mathrm{2}}+{\beta }^{\mathrm{2}}{p}_{{\tau }_{\mathrm{2}}}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}+{\mu }_{\mathrm{1}}{p}_n=\mathrm{0}⟹\left|\mathrm{\Delta }{u}_{\tau }\right|\ge \mathrm{0};\end{array}$

(5)

if

$\begin{array}{c}\left|{\mathbf{p}}_{\tau }\right|<-{\mu }_f{p}_n⟹\left|d{u}_{\tau }^{\left(\mathrm{1}\right)}-d{u}_{\tau }^{\left(\mathrm{2}\right)}\right|=\mathrm{0},\\ \left|{\mathbf{p}}_{\tau }\right|=-{\mu }_f{p}_n⟹\left|d{u}_{\tau }^{\left(\mathrm{1}\right)}-d{u}_{\tau }^{\left(\mathrm{2}\right)}\right|\ge \mathrm{0};\end{array}\}\mathrm{on}\hspace{0.17em}\hspace{0.17em}{S}_c,$

(6)

where all the variables and parameters are notated conventionally.

2.2. Constitutive Equations and Orthotropic Friction Law

The constitutive equation can be obtained by applying Taylor expansion to (3):

$f_{\alpha }^{\mathrm{0}}+{\mathbf{W}}_{\alpha }d\epsilon -M_{\alpha }\lambda \le \mathrm{0},$

(7)

where α = 1, 2,…, y, λ ≥ 0, f_α⁰ is the initial value of the yield function at instant, and

$\begin{array}{c}{\mathbf{W}}_{\alpha }=\left(\frac{\partial f_{\alpha }}{\partial \mathit{\sigma }}\right)\mathbf{D},\\ M_{\alpha }={\mathbf{W}}_{\alpha }{\left(\frac{\partial g}{\partial \mathit{\sigma }}\right)}^T-{\left(\frac{\partial f_{\alpha }}{\partial {\mathbf{\epsilon }}^p}\right)}^T{\left(\frac{\partial g}{\partial \mathit{\sigma }}\right)}^T+\frac{\partial f_{\alpha }}{\partial {\mathbf{\kappa }}_{\alpha }}{\mathbf{h}}^T.\end{array}$

(8)

The boundary condition on the elastic-plastic contact surface S _c (Figure 2) can also be written otherwise. Particularly, the Coulomb frictional law for 3D contact problems that is not expressed identically to that for 2D problems can be expressed as

$\begin{array}{c}{C}_{\infty }\left(p_n\right)=\left\{{\mathbf{p}}_{\tau }:f=\left|{\mathbf{p}}_{\tau }\right|+{\mu }_f{p}_n\le \mathrm{0}\right\},\\ {\mathbf{p}}_{\tau }=p_{{\tau }_{\mathrm{1}}}\mathbf{i}+p_{{\tau }_{\mathrm{2}}}\mathbf{j}.\hspace{0.17em}\end{array}$

(9)

Figure 2:

Two contact bodies.

To establish the complementary equation, (9) can be linearized as

$\begin{array}{c}{C}_{N_f}=\left\{{\mathbf{p}}_{\tau }:f_i\left({\mathbf{p}}_{\tau },p_n\right)\le \mathrm{0},\hspace{0.17em}i=\mathrm{1,2},\dots ,N_f\right\}\\ f_i=\left[\cos {\alpha }_i,\sin {\alpha }_i,{\mu }_f\right]{\left[p_{{\tau }_{\mathrm{1}}},p_{{\tau }_{\mathrm{2}}},p_n\right]}^T.\end{array}$

(10)

Introducing penalty factors E_τ, E _n in the tangential and normal directions of the contact surface and decomposing the contact strains into an elastic part (describing microslip) and a plastic part (describing macroslip), the following equations can be obtained:

$d{\mathbf{\epsilon }}_c=d{\mathbf{\epsilon }}_c^e+d{\mathbf{\epsilon }}_c^p,$

(11)

$d{\mathbf{p}}_c={\mathbf{D}}_c\left(d{\mathbf{\epsilon }}_c-d{\mathbf{\epsilon }}_c^p\right),$

(12)

where

$\begin{array}{c}d{\mathbf{\epsilon }}_c={\left[d{\epsilon }_{{\tau }_{\mathrm{1}}},d{\epsilon }_{{\tau }_{\mathrm{2}}},d{\epsilon }_n\right]}^T,\\ d{\mathbf{p}}_c={\left[d{p}_{{\tau }_{\mathrm{1}}},d{p}_{{\tau }_{\mathrm{2}}},d{p}_n\right]}^T,\\ {\mathbf{D}}_c=\left[\begin{array}{ccc}{E}_{\tau }& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& E_{\tau }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& E_n\end{array}\right],\\ d{\epsilon }_{{\tau }_{\mathrm{1}}}=d{u}_{{\tau }_{\mathrm{1}}}^{\left(\mathrm{1}\right)}-d{u}_{{\tau }_{\mathrm{1}}}^{\left(\mathrm{2}\right)}=\mathrm{\Delta }{u}_{{\tau }_{\mathrm{1}}},\\ d{\epsilon }_{{\tau }_{\mathrm{2}}}=d{u}_{{\tau }_{\mathrm{2}}}^{\left(\mathrm{1}\right)}-d{u}_{{\tau }_{\mathrm{2}}}^{\left(\mathrm{2}\right)}=\mathrm{\Delta }{u}_{{\tau }_{\mathrm{2}}},\\ d{\epsilon }_n=d{u}_n^{\left(\mathrm{1}\right)}-d{u}_n^{\left(\mathrm{2}\right)}+{\delta }^{*}=\mathrm{\Delta }{u}_n+{\delta }^{*},\end{array}$

(13)

and δ^* is the initial normal gap between contact surfaces.

In analogy with the expression of elastic-perfectly plastic constitutive model, slip function ${\tilde{f}}_i$ and slip potential function ${\tilde{g}}_i$ , which are given as follows, are introduced to describe the contact:

$\begin{array}{c}{\tilde{f}}_i=p_{{\tau }_{\mathrm{1}}}\cos {\alpha }_i+p_{{\tau }_{\mathrm{2}}}\sin {\alpha }_i+{\mu }_f{p}_n\le \mathrm{0},\\ {\tilde{g}}_i=p_{{\tau }_1}\cos {\alpha }_i+p_{{\tau }_2}\sin {\alpha }_i-c,\\ {\tilde{f}}_{N_f+1}=p_n\le \mathrm{0},\\ {\tilde{g}}_{N_f+1}=p_n,\end{array}$

(14)

where i = 1, 2,…, N _f , c is a constant and the terms with “∼” are the variables describing contact. Therefore, the contact plastic (slip or separation) incremental displacement vector can be acquired by

$d{\mathbf{\epsilon }}_c^p={{{\displaystyle \sum }}_{i=\mathrm{1}}^{N_f+\mathrm{1}}{\tilde{\lambda }}_i\frac{\partial {\tilde{\mathbf{g}}}_i}{\partial {\mathbf{p}}_c}={\left(\frac{\partial \tilde{\mathbf{g}}}{\partial {\mathbf{p}}_c}\right)}^T}^{}\tilde{\mathbf{\lambda }},\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}$

(15)

where $\tilde{g}={[{\tilde{g}}_{\mathrm{1}},{\tilde{g}}_{\mathrm{2}},\dots ,{\tilde{g}}_{N_f+\mathrm{1}}]}^T$ and $\tilde{\lambda }={[{\tilde{\lambda }}_{\mathrm{1}},{\tilde{\lambda }}_{\mathrm{2}},\dots ,{\tilde{\lambda }}_{N_f+\mathrm{1}}]}^T\ge \mathrm{0}$ . Thus, the contact boundary conditions are changed into

${\tilde{f}}_i^0+{\tilde{\mathbf{W}}}_{i}d{\epsilon }_c-{\tilde{\mathbf{M}}}_i\tilde{\mathbf{\lambda }}+{\tilde{v}}_i=\mathrm{0},$

(16)

where

$\begin{array}{c}{\tilde{\lambda }}_i\ge \mathrm{0},\hspace{1em}\hspace{1em}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\tilde{v}}_i\ge \mathrm{0},\hspace{1em}i=\mathrm{1,2},\dots ,N_f+\mathrm{1},\\ {\tilde{\mathbf{W}}}_i={\left[\frac{\partial {\tilde{f}}_i}{\partial {\mathbf{p}}_c}\right]}^T{\mathbf{D}}_c\begin{array}{cc},\hspace{1em}\hspace{1em}& \tilde{\mathbf{M}}\end{array}={\left[\frac{\partial {\tilde{f}}_i}{\partial {\mathbf{p}}_c}\right]}^T{\mathbf{D}}_c{\left(\frac{\partial \tilde{\mathbf{g}}}{\partial {\mathbf{p}}_c}\right)}^T.\end{array}$

(17)

2.3. PQP Solution

Among the admissible displacement fields, the actual solution minimizes the potential energy

$\begin{array}{l}\text{Π}\left[\mathit{\lambda },\tilde{\mathit{\lambda }}\right]\\ \hspace{1em}={{\displaystyle \int }}_{\mathrm{\Omega }}^{}\left[\frac{d{u}_{i,j}{D}_{ijkl}d{u}_{k,l}}{2}-{\lambda }_{\alpha }{R}_{kl\alpha }d{u}_{k,l}-d{b}_{i}d{u}_i\right]d\mathrm{\Omega }\\ \hspace{0.17em}\hspace{1em}+{{\displaystyle \int }}_{s_c}^{}\left(\frac{d{\epsilon }_{ci}{D}_{cij}d{\epsilon }_{ij}}{2}-{\tilde{\lambda }}_{\alpha }{\tilde{R}}_{k\alpha }d{\epsilon }_{ck}\right)ds-{{\displaystyle \int }}_{s_p}^{}d{\stackrel{-}{p}}_{i}d{u}_{i}ds,\end{array}$

(18)

controlled by (7) and (16), where

$\begin{array}{c}{R}_{kl\alpha }=\left(\frac{\partial g_{\alpha }}{\partial {\sigma }_{ij}}\right)D_{ijkl},\\ {\tilde{R}}_{k\alpha }=\left(\frac{\partial {\tilde{g}}_{\alpha }}{\partial p_{ci}}\right)D_{cik},\\ {\epsilon }_c=\mathrm{\Delta }{u}_c,\\ D_{cij}={\delta }_{ij}{E}_{\tau }={\delta }_{ij}{E}_n={\delta }_{ij}E.\end{array}$

(19)

Posterior to introducing the interpolation shape function matrix N and assembling the total energy, (12) can be derived as

$\text{Π}\left[{\mathbf{\lambda }}^{\prime }\right]=\frac{\mathrm{1}}{\mathrm{2}}d{\mathbf{u}}^T\mathbf{K}d\mathbf{u}-d{\mathbf{u}}^T\left(\mathrm{\Phi }{\mathbf{\lambda }}^{\prime }+{\mathbf{p}}_{\mathrm{0}}-{\mathbf{p}}_{\delta }\right),$

(20)

where

$\begin{array}{l}\mathbf{K}={{\displaystyle \int }}_{\mathrm{\Omega }}^{}{\mathbf{B}}^T\mathbf{D}\mathbf{B}\hspace{0.17em}d\mathrm{\Omega }+{{\displaystyle \int }}_{s_c^{}}^{}{\mathbf{N}}^T{\mathbf{D}}_c\mathbf{N}\mathbf{\hspace{0.17em}}ds,\\ {\mathbf{p}}_{\mathrm{0}}={{\displaystyle \int }}_{\mathrm{\Omega }}^{}{\mathbf{N}}^{T}d\mathbf{b}\mathbf{\hspace{0.17em}}d\mathrm{\Omega }+{{\displaystyle \int }}_{s_p^{}}^{}{\mathbf{N}}^{T}d\stackrel{-}{\mathbf{p}}\hspace{0.17em}ds,\\ {\mathbf{p}}_{\delta }={{\displaystyle \int }}_{s_c^{}}^{}{\mathbf{N}}^T{\mathbf{D}}_c\mathbf{N}\mathbf{\hspace{0.17em}}ds{\mathit{\delta }}_c^{*},\hspace{0.17em}\hspace{0.17em}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\mathbf{\lambda }}^{\prime }={\left[{\mathbf{\lambda }}^T,{\tilde{\mathbf{\lambda }}}^T\right]}^T,\\ \mathrm{\Phi }={{\displaystyle \int }}_{\mathrm{\Omega }}^{}{\mathbf{B}}^T{\mathbf{R}}^{T}d\mathrm{\Omega }+{{\displaystyle \int }}_{s_c^{}}^{}{\mathbf{N}}^T{\tilde{\mathbf{R}}}^{T}ds.\end{array}$

(21)

The control (7) and (16) can be derived as

$\begin{array}{l}\mathbf{C}\mathbf{\hspace{0.17em}}d\mathbf{u}-\mathbf{U}{\mathbf{\lambda }}^{\prime }-({\mathbf{t}}_{\mathrm{0}}+{\mathbf{t}}_{\delta })+{\mathbf{v}}^{\prime }=\mathrm{0},\\ {{\mathbf{v}}^{\prime }}^T{\mathbf{\lambda }}^{\prime }=\mathrm{0},\end{array}$

(22)

where

${\mathbf{v}}^{\prime },{\mathbf{\lambda }}^{\prime }\ge \mathrm{0},$

(23)

$\begin{array}{l}\mathbf{C}={{\displaystyle \int }}_{\mathrm{\Omega }}^{}\mathbf{W}\mathbf{B}\hspace{0.17em}d\mathrm{\Omega }+{{\displaystyle \int }}_{{\text{s}}_{\text{c}}^{}}^{}\tilde{\mathbf{W}}\mathbf{N}\hspace{0.17em}ds,\\ \mathbf{U}={{\displaystyle \int }}_{\mathrm{\Omega }}^{}\mathbf{M}\hspace{0.17em}d\mathrm{\Omega }+{{\displaystyle \int }}_{s_c^{}}^{}\tilde{\mathbf{M}}\hspace{0.17em}ds,\\ {\mathbf{t}}_{\delta }=-{{\displaystyle \int }}_{s_c^{}}^{}{\tilde{\mathbf{W}}}_{e}ds{\mathit{\delta }}_c^{*},\\ {\mathbf{t}}_{\mathrm{0}}=-{{\displaystyle \int }}_{\mathrm{\Omega }}^{}{\mathbf{f}}_{\mathrm{0}}^{}d\mathrm{\Omega }-{{\displaystyle \int }}_{s_c^{}}^{}{\tilde{\mathbf{f}}}_{\mathrm{0}}^{}ds.\end{array}$

(24)

From (20) and (22), the following equations can be obtained:

$\begin{array}{l}{\mathbf{v}}^{\prime }-(\mathbf{U}-\mathbf{C}{\mathbf{K}}^{-\mathrm{1}}\mathrm{\Phi }){\mathbf{\lambda }}^{\prime }\\ \hspace{1em}=-\mathbf{C}{\mathbf{K}}^{-\mathrm{1}}{\mathbf{p}}_{\mathrm{0}}+{\mathbf{t}}_{\mathrm{0}}-(\mathbf{U}-\mathbf{C}{\mathbf{K}}^{-\mathrm{1}}\mathrm{\Phi }){\mathit{\delta }}_c^{*},\\ {{\mathbf{v}}^{\prime }}^T{\mathbf{\lambda }}^{\prime }=\mathrm{0},\end{array}$

(25)

where v_′,λ_′ ≥ 0, and δ _c ^* is expanded to the space of λ_′.

It is necessary to eliminate the penalty factor E in the standard linear complementarity problem (25) to correct the solution. The matrices in (25) are in the forms of

$\begin{array}{l}\mathbf{K}=\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{11}}& {\mathbf{K}}_{\mathrm{12}}\\ {\mathbf{K}}_{\mathrm{21}}& {\mathbf{K}}_{\mathrm{22}}+E{{{\mathbf{K}}_{\mathrm{22}}^{\prime }}^{}}_{}\end{array}\right],\\ {\mathbf{p}}_{\mathrm{0}}=\left[\begin{array}{c}{\mathbf{p}}_{0i}\\ {\mathbf{p}}_{0o}\end{array}\right],\hspace{1em}\hspace{1em}{\mathbf{t}}_{\mathrm{0}}=\left[\begin{array}{c}{\mathbf{t}}_{0i}\\ {\mathbf{t}}_{0o}\end{array}\right],\\ \mathbf{C}=\left[\begin{array}{cc}{\mathbf{C}}_{\mathrm{11}}& {\mathbf{C}}_{\mathrm{12}}\\ \mathrm{0}& E{{\mathbf{C}}_{\mathrm{22}}^{\prime }{}^{}}_{}\end{array}\right],\hspace{1em}\hspace{1em}\mathbf{U}=\left[\begin{array}{cc}{\mathbf{U}}_{\mathrm{11}}& \mathrm{0}\\ \mathrm{0}& \mathrm{\hspace{0.17em}\hspace{0.17em}}E{{{\mathbf{U}}_{\mathrm{22}}^{\prime }}^{}}_{}\end{array}\right],\\ \begin{array}{c}𝚽=\left[\begin{array}{cc}{𝚽}_{\mathrm{11}}& \mathrm{0}\\ {𝚽}_{\mathrm{21}}& \mathrm{\hspace{0.17em}\hspace{0.17em}}E{{{𝚽}_{\mathrm{22}}^{\prime }}^{}}_{}\end{array}\right].\mathrm{\hspace{0.17em}\hspace{0.17em}}\end{array}\end{array}$

(26)

After a series of derivations, (25) can be divided into

$\begin{array}{c}\mathbf{v}-({\mathbf{U}}_{\mathrm{11}}-{\mathbf{C}}_{\mathrm{11}}{\mathbf{K}}_{\mathrm{11}}^{-\mathrm{1}}{𝚽}_{\mathrm{11}})\mathbf{\lambda }-({\mathbf{C}}_{\mathrm{11}}{\mathbf{K}}_{\mathrm{11}}^{-\mathrm{1}}{\mathbf{K}}_{\mathrm{12}}-{\mathbf{C}}_{\mathrm{12}}){{{𝚽}_{\mathrm{22}}^{\prime }}^{}}_{}\tilde{\mathbf{\lambda }}\\ \hspace{1em}\hspace{1em}=-{\mathbf{C}}_{\mathrm{11}}{\mathbf{K}}_{\mathrm{11}}^{-\mathrm{1}}{\mathbf{p}}_{0i}+{\mathbf{t}}_{0i}-({\mathbf{C}}_{\mathrm{11}}{\mathbf{K}}_{\mathrm{11}}^{-\mathrm{1}}{\mathbf{K}}_{\mathrm{12}}-{\mathbf{C}}_{\mathrm{12}}){\mathit{\delta }}_c^{*},\\ \tilde{\mathbf{v}}-{{{\mathbf{C}}_{\mathrm{22}}^{\prime }}^{}}_{}({\mathbf{K}}_{\mathrm{21}}{\mathbf{K}}_{\mathrm{11}}^{-\mathrm{1}}{𝚽}_{\mathrm{11}}-{𝚽}_{\mathrm{21}})\mathbf{\lambda }-{{{\mathbf{C}}_{\mathrm{22}}^{\prime }}^{}}_{}{\mathbf{K}}_{\mathrm{22}}^{*}{{{𝚽}_{\mathrm{22}}^{\prime }}^{}}_{}\tilde{\mathbf{\lambda }}\\ \hspace{1em}\hspace{1em}=-{\mathbf{C}}_{\mathrm{22}}^{\prime }{\mathbf{p}}_{\mathrm{0}}^{*}+{\mathbf{t}}_{0o}-{\mathbf{C}}_{\mathrm{22}}^{\prime }{\mathbf{K}}_{\mathrm{22}}^{*}{\mathit{\delta }}_c^{*},\\ {\mathbf{v}}^T\mathbf{\lambda }=\mathrm{0},\hspace{1em}\hspace{1em}{\tilde{\mathbf{v}}}^T\tilde{\mathbf{\lambda }}=\mathrm{0},\end{array}$

(27)

where v, λ, $\tilde{\mathbf{v}}$ , $\tilde{\mathbf{\lambda }}\ge \mathrm{0}$ ,

${\mathbf{K}}_{\mathrm{22}}^{*}={\mathbf{K}}_{\mathrm{22}}-{\mathbf{K}}_{\mathrm{21}}{\mathbf{K}}_{\mathrm{11}}^{-\mathrm{1}}{\mathbf{K}}_{\mathrm{12}},\hspace{1em}\hspace{1em}{\mathbf{p}}_{\mathrm{0}}^{*}={\mathbf{p}}_{0o}-{\mathbf{K}}_{\mathrm{21}}{\mathbf{K}}_{\mathrm{11}}^{-\mathrm{1}}{\mathbf{p}}_{0i}.$

(28)

Apparently, the expressions in (28) duplicate the load vector and the stiffness matrix in substructuring which can thus be used to deal with contact problems, and the algorithm herein can be readily implemented in a code.

3. Details of the Computational Model

The compressor impeller, which comprises identical segments circumferentially, gives rise to a rotationally periodic structure. The impeller, with a shaft hole of diameter 34 mm, is mounted onto compressor shaft and sleeve via interference fit. The material of the sleeve and axle 7 is steel and that of the impeller is aluminum. Although adjacent blades are not largely gapped, impeller requires to be computed by using the three-dimensional model other than an axial symmetric one because of the cyclic symmetric structure that should be analyzed by substructuring. The global stiffness matrix is created merely by dividing meshes into one substructure mode. Besides, the computed results can be applied to other similarly constructed superelements, which thus significantly lower the computation cost. In other words, increasing repeated structures will augment the corresponding computational efficiency.

Table 1 summarizes the overall dimensions and the basic design data of the impeller, and the physical characteristics of LD7 are listed in Table 2.

Table 1:

Major dimensions of the impeller.

Basic design data	Parameters
Impellertip diameter d1/mm	310
Inducertip diameter d2/mm	175
Ratedspeed w/rpm	25000
Shaftdiameter d3/mm	34
Shaftsleeve width h/mm	7
Numberof blades z	24
Material	Aluminumalloy (LD7)

Table 2:

Physical characteristic data of LD7.

Density (kg·m−3)	2800	σ b (Mpa)	338
Extensibility	10	σ S (Mpa)	261
Poisson ratio	0.33	Elastic modulus (Gpa)	71
		shear modulus (Gpa)	2.68

The substructure model is exhibited in Figure 3. The computational models of the impeller are shown in Figure 4.

Figure 3:

The model of substructure.

Figure 4:

The whole computational model of impeller.

The computational model uses a four-node tetrahedron solid element with three degrees of freedom at each node to truly reflect the impeller. The influences of temperature and pressure loading are neglected because they are smaller than the centrifugal force due to high-speed rotation. Meanwhile, the centrifugal stiffening effect [18] is also considered. In the FE contact analysis, contact surfaces are defined to allow nodal sliding.

Simulating the tested constitutive curve reveals that the elastoplastic computation of the impeller follows a linear hardening model (Figure 5).

Figure 5:

Stress-strain tested curve of LD7.

To evaluate the changes of displacements and stresses regarding the impeller, 60% of the full load is first applied according to the results of the elastic model with the maximum induced stress close to the proportional limit, after which the impeller undergoes 60, 10, 10, 10, and 10% of the full load, unloading, and then reloading at 25000 rpm.

The compressor essentially succumbs to centrifugal forces owing to the high-speed rotation in the midst of overspeed and operation. δ₁ and δ₂ are defined as the amount of interference between shaft sleeve and impeller as well as that between shaft and shaft sleeve, respectively. As to δ₁ (0.03, 0.04, and 0.05 mm) and δ₂ (0.04, 0.05, and 0.06 mm), the load cases for different wall thicknesses h (7, 6, 5, and 4 mm) and overspeed rpms w (30000, 31000, 32000, 33000, 34000, and 35000 rpm) are computed.

4. Results and Discussions

4.1. Influence of Overspeed

The solid curves in Figure 6 represent the distribution of contact stresses of the external surface of shaft sleeve (surface between shaft sleeve and impeller), whereas the dashed curves correspond to that of the internal surface (surface between shaft and shaft sleeve). Their distributions along the axial direction at five rotational speeds are also illustrated (δ₁ = 0.03 mm, δ₂ = 0.04 mm). The right edge of the shaft sleeve external surface is of higher contact stress than other regions because impeller and shaft sleeve do not keep contacting persistently. The impeller is suspended at x = 47.0 mm (Figure 3) without contacting with the shaft sleeve external surface, the root discontinuity geometry which generates a stress concentration. Encountering the same amount of interference, the contact stresses of contacting surface reduce with accelerating rotation because of the raised centrifugal forces. In addition, the contact stresses of the surface edge barely change, and those of internal and external surfaces change identically.

Figure 6:

(a) Normal contact stress distributions of the external surface of the shaft sleeve along the axial direction under different overspeed rpms. (b) Normal contact stress distributions of the internal surface of the shaft sleeve along the axial direction under different overspeed rpms.

Half of the contact points between impeller and shaft sleeve are separated at the rotational speed of 30000 rpm (Figure 6).

The solid curves in Figure 7 represent the distribution of contact stresses of the shaft sleeve external surface. The normal displacements of contact point pairs between shaft sleeve and impeller increase along the axial direction with rising radial dimension of the impeller (Figure 7).

Figure 7:

Variation of normal displacement under different overspeed rpms.

The normal displacement values for x = − 30.25∼ − 43.25 mm at the contact surface between shaft sleeve and impeller peak and thereafter decrease down to 0 along the axial direction. The contact stresses of contacting surfaces are inevitably nonuniformly distributed owing to the nonuniform distribution of normal displacements. We postulate that the nonuniform initial amount of interference manages to avoid displacement and unify contact stress. For instance, the amount of interference linearly rises (declines). Figure 7 also suggests that the normal displacement of contact point pairs is elevated with accelerating rotation.

4.2. Influence of Amount of Interference

The solid curves in Figure 8 correspond to the distribution of contact stresses of the shaft sleeve external surface, while the dashed ones represent that of the internal surface. The relationships between contact stresses and rotational speeds for various amounts of interference (case 1: δ₁ = 0.03, δ₂ = 0.04; case 2: δ₁ = 0.04, δ₂ = 0.05; case 3: δ₁ = 0.05, δ₂ = 0.06) are derived at the central contact point x = 21.25 mm. Figure 8 shows that the contact stresses of shaft sleeve external and internal surfaces linearly drop with accelerating rotation, and they also decrease with reducing amount of interference at the same rotational speed. The contact stress of internal surface which exceeds that of external one also decreases more rapidly, which can be attributed to the difference between materials (shaft sleeve and shaft: steel; impeller: aluminum) and the remarkably lower rigidity of aluminum than that of steel.

Figure 8:

Variation of normal contact stress with overspeed rpm when the amount of interferences is different.

Consequently, the contact stresses of impeller, shaft, and shaft sleeve are essentially affected by the amount of interference between shaft sleeve and shaft which prevents the contact stress there between from fast decrease at high rotational speed, thereby guaranteeing the safety of overspeed prestressing.

4.3. Influence of the Wall Thickness of Shaft Sleeve

The solid curves in Figure 10 represent the distribution of contact stresses of the shaft sleeve external surface, while the dashed ones correspond to that of the internal surface. The relationships between contact stresses and rotational speeds for various wall thicknesses of shaft sleeve at the central contact point x = 21.25 mm are plotted in Figure 9. The contact stresses of both external and internal surfaces are lowered with accelerating rotation. With increasing wall thickness, however, the contact stress of internal surface increases, while that of external surface changes conversely at the same rotational speed. Thickening the wall of shaft sleeve that depends on the radius of external surface indicates increased impeller internal diameter and decreased impeller stiffness in case of constant amount of interference between shaft and shaft sleeve. Thus, the external surface is of lower contact stress, and we can ascribe the elevated contact stress of internal surface to the same reason.

Figure 9:

Variation of normal contact stress with overspeed when the wall thicknesses of shaft sleeve are different.

Figure 10:

Variation of maximum residual Von Mises stress of impeller under different overspeed rpms.

4.4. Maximum Residual Stress and Displacement

As to the impeller, Figures 10 and 11 exhibit the dependences of maximum residual Von Mises stress and displacement on overspeed rpm. The maximum residual Von Mises stress at 30000 rpm∼35000 rpm is located close to the internal hole rear of the compressor. The maximum residual radial displacement resides near the inner hole rear of the impeller, while the axial one is located adjacent to the disk edge. The maximum residual Von Mises stress and displacement both rise rapidly with accelerating rotation, which accords with the continuous expansion of plastic rear zone.

Figure 11:

Variation of maximum residual displacement of impeller under different overspeed rpms.

4.5. Stress and Displacement under Working Conditions

Figure 12 presents the dependence of the maximum Von Mises stress of the impeller on overspeed rpm, in which the maximum Von Mises stress at working speed declines slowly with rising rotational speeds after overspeed prestressing, suggesting that the impact of prestressing on the impeller is weakened with accelerating rotation.

Figure 12:

Variation of maximum Von Mises stress of impeller under different overspeed rpms.

Figure 13 exhibits the dependence of the maximum displacements of the impeller on overspeed rpm. The axial and radial displacements, which increase with elevating rotational speeds overall, first increase moderately at the speeds lower than 33000 rpm and then rapidly at those exceeding 33000 rpm.

Figure 13:

Variation of maximum displacement of impeller under different overspeed rpms.

5. Conclusions

1. The PQP protocol manages to solve convergence problems and to reduce sophisticated high nonlinear problems (e.g., 3D elastoplastic problems) to standard linear complementarity ones by approximating the friction law via a polyhedral one.

2. At the overspeed of 30000 rpm, the contact region of the impeller radically sized higher is susceptible to displacement, and half of the contact points between shaft sleeve and impeller may be gapped. The computation results show that at different overspeed rpms, the impeller mass distributions and the normal relative displacements of the contact point pairs between shaft sleeve and impeller are apparently associated.

3. The contact stresses of shaft-shaft sleeve-impeller are mainly affected by the amount of interference between shaft and shaft sleeve. The contact stresses of the internal and external surfaces of shaft sleeve are linearly lowered with accelerating rotation, and the internal surface is of higher contact stress than the external one. To qualify the overspeed prestressing, it is critical to strictly control the amount of interference between shaft and shaft sleeve to prevent the contact stress from fast increasing. Meanwhile, the nonuniform initial amount of interference can get rid of displacement and unify contact stress.

4. The contact stress of the internal surface of shaft sleeve rises and that of the external surface reduces with thickening wall. The wall thickness should be minimized to maintain sufficient contact stress and to prevent it from fast changes along the axial direction.

5. The effects of prestressing minimizing the maximum impeller stress at working speed should be considered in identifying the optimal overspeed rpm. In the meantime, the impeller succumbs to tolerable residual deformation after reasonably costing overspeed. In this study, the impeller (calculation example) functions optimally at the overspeed rpms from 31000 rpm to 33000 rpm at which the maximum stress decreases to 167 Mpa, and the maximum displacement at working speed is rapidly augmented at over 33000 rpm.