© J. Li et al., published by EDP Sciences 2022

1 Introduction

Conventional steam Rankine cycle operate at a temperature of less than 550 °C [1–3], based on a trade-off between material selection, safety and reliability of circulatory mediator and rotor speed of the turbomachinery accounting also for bearings and seals induced restrictions of operating temperatures and pressures [4–6]. In recent years, reliability combined with the high thermal efficiencies related advantages has led to an increasing attention to supercritical carbon dioxide (SCO₂) energy conversion system [7]. As a potential alternative power conversion system, The SCO₂ Brayton cycle has been extensively studied [8]. The relevance of inlet temperature is reflected by the numerous numerical and experimental studies that have been implemented on this subject in the past few years. To meet the demand for high thermal efficiency, it is widely acknowledged that heat sources ranging between 500 °C and 900 °C and supercritical (>30.98 °C and >7.38 MPa) conditions will be required. Therefore, elevated temperature turbines may still remain a pivotal component within future SCO₂ power cycles. On the other hand, in the drive towards compact physical footprint, and to surmount challenges of the rotor dynamics related to high speed, attention has turned to the small scale of thermodynamic turbine [9]. As the working speeds (300 kW, 50,000 rpm) and inlet temperature of the small SCO₂ turbines are often extreme, a strict optimization design of the shaft bearing system is required. The integrated optimization for multiple physics including thermomechanics, rotor-bearing dynamics.

Researches on these optimization algorithms for bearing-rotor system applications are neither unique nor trivial as they have applications in the rotating machinery within turbine engines and synchronous motors. The development of contemporary designs of bearing-rotor systems depends on the amelioration of methods of simulation and design as well as algorithms of optimization. Combined with finite element analysis (FEA), it allows the time-consuming repetitive work of optimization of prototypes to be avoided. Earlier, Hamidreza et al. [10] proposed a design method of flexible rotor support parameters to minimize the force transmission in the vicinity of the first critical speed in a preliminary design phase. Detailed optimum variables were used to design the support flexibility and damping. However, this approach shows limitations with regard to the optimization of multi-objective functions due to the deficiencies of the H _∞ and H ₂ algorithms. SangJoon et al. [11] conducted an optimization study on the rotor structure of a compound rotorcraft. The objective functions are rotor weight and vibration performance. Furthermore, comprehensive optimization analysis were conducted by Omar [12] using rigorous multi-objective optimization algorithm in the electric machinery model to make a trade-off between rotor geometries and loss objectives.

Over the past decades, natural heuristic algorithms based on statistics have been extensively developed [13]. Genetic algorithm (GA) is one of the most commonly used. Angie et al. [14] carried out a CFD design optimization of the radial-inflow turbine with respect to efficiency and power output improvement. For this purpose, they used Gaussian radial basis function in combination with an improved genetic algorithm (NSGA II) to optimize the parameterized shape of the turbine rotor system. Hybrid genetic algorithm investigations performed by Copeland et al. [15] sought to couple an evolutionary algorithm with a surrogate model in the application of aerodynamic shape optimization aiming at enhance micro-wave rotor turbine efficiency. In addition to wave rotors, a hybrid of genetic algorithm with the Nelder-Mead method have been used in the preliminary design of organic rankine cycle (ORC) axial turbine to improve the flow efficiency [16]. Another use of genetic algorithm was presented in Joseph's [17] investigation, which coupled with the controlled elitist multi-objective optimization method to evaluate the unbalance response performances of small turbine engine rotor system. The results of comparing the amplitude of vibrations show that the optimum design model is superior to constant Rayleigh damping model.

Many investigations are conducted focusing on torque and ripple of motor rotor from the perspective of topological structure by using hybrid optimization algorithm (HOA). A typical HOA is a mixed mathematical statistical method of topology optimization (TO) and genetic algorithm (GA), which provides correlations between the non-parameterize shape and the design variables. Shingo et al. [18] applied a hybrid algorithm to optimize the rotor shape and topology of the flux barrier. A feasible trade-off between the rotor torque and ripple is proposed with respect to pivotal parameters that affect the problem. Nag and Harish et al. [19] investigated interference parameters of an Orbit motor rotor system. Their main objectives were to determine interference effects between the rotor and stator to identify further unbalanced torque reduction potential by the application of GA. Regarding topology optimization, Takeo et al. [20] conducted numerical investigations aimed at optimizing the rotor topology of a PM synchronous motors by using genetic algorithms. The optimized topology showed a 30% increase in torque compared to the values obtained by the experimental motor. Karol et al. [13] also conducted a similar optimization study using the gray wolf metaheuristic algorithm, which considered the functional parameters of the motor related to energy efficiency.

The probabilistic methods are particularly effective in optimizing a single objective function [21]. By improving the search efficiency, increasingly attention-grabbing bionic stochastic algorithms, can also iterate to the extremums of an objective function. Particle swarm optimization (PSO) is one of the most widely known swarm intelligence bionic algorithms. Recent research studies have taken up territory particle swarm optimization (TPSO) algorithm and applied it to the optimization design of permanent magnet synchronous motor (PMSM) to reduce cogging torque [22]. Mutra [23] investigated interference effects between the rotor spin speed and the foil bearing parameters (including bearing clearance, bump foil thickness, and foil pitch) of a conventional high-speed turbocharger rotor system by supervised neural network scheme in order to identify the nonlinear relationship between the bearing parameters and vibration response data. Furthermore, an optimal design of foil bearing parameters is carried out using an improved particle swarm optimization algorithm. In the damage identification method, the bat optimization algorithm (BA) combined with multiaxial vibration data investigations were performed for a helicopter main rotor inverse damage location [24]. One study conducted by Karol et al. used the gray wolf algorithm (GWA) to optimize variables related to the motor rotor structure to design a low-power line-start synchronous motor [13]. The results show that up to 1.53% efficiency increase was obtained, in comparison to the induction motor.

In the pre-design of the bearing rotor system, numerical calculation is usually used to optimize the relevant parameters (such as the stiffness, span and damping of the support). Finite element analysis (FEA) provides the engineer with a numerical implementation and verification of the prototypes. Considering an advanced SCO₂ turbine design, there are two feasible ways to improve the performance of the turbine rotor, namely to increase the inlet temperature and speed. Therefore, the attenuation of harmonic response amplitude and thermal management has always been an important topic during the SCO₂ rotor-bearing systems development. The present work is part of the SCO₂ turbine system design optimization project, which is well aligned with the required innovation for high-performance turbine. The project is aimed on the rotor dynamics design optimization and vibration control of a small scale SCO₂ turbine rotor system. Most of these previous investigations utilized single-objective and nature-inspired algorithms to solve the parameters and topology of rotor-bearing system. Nature-inspired algorithms (non-deterministic algorithms) are significantly affected by subjective factors and make parametric studies cumbersome. Moreover, multiple Pareto solutions of multivariate optimization are difficult to capture, especially in the region of multiple conflicting objectives. For this purpose, a multi-objective and multi-scale quantum harmonic oscillator algorithm (mMQHOA) combined with transfer matrix method is developed based on standard MQHOA [25]. Our investigations improved the previous algorithm to consider more than two kinds of objectives, namely, resonance amplitude, critical speed and its position, and designed the rotor structure considering the cooling system.

The bearing-rotor system of a small scale SCO₂ turbine, which is shown in Figure 1, was rotor-dynamically optimized within the present work. Initially, this paper presents the optimization of the rotor support system considering a journal cooler as well as two types of bearings (including angular contact ball bearings (ACBB) and Oil-lubricated tilted pad bearings (OTPB)). Bearings and cantilever impeller are a major source of vibration in SCO₂ turbine rotor systems, producing a demand for strict vibration evaluation and frequency response analysis. The application of mMQHOA in the dynamic design is to minimize the resonance response amplitude and adjust the position of the critical frequency relative to the operating frequency. In order to determine the advantages of the optimized rotor system, it is compared to a reference structure, which was supported by four-pad oil film bearing during the preliminary design stage of the FE demonstrator. Then, the finite element method is further used to evaluate the sensitivity of the rotor system's harmonic response to parameters such as bearing span and cooler film stiffness.

Fig. 1 A schematic view of the investigated bearing-rotor system (a) and its cooler (b). (a) Schematic diagram of bearing-rotor system. (b) Schematic diagramc of the cooler.

2 Methodology

2.1 Structure of the rotor-bearing system

The force coefficients of the bearing are combined with the mass, damping and stiffness of the hub to form a composite matrix, as shown below$$K_e=\frac{\left[K_{oil}+A{C}_{oil}\right]\left[K_{hub}+A{C}_{\mathrm{hub}}+A^2{M}_{\mathrm{hub}}\right]}{\left[K_{oil}+A{C}_{oil}+K_{\mathrm{hub}}+A{C}_{\mathrm{hub}}+A^2{M}_{\mathrm{hub}}\right]}$$(1)

This is the dynamic model of a series spring-mass-damping system. The K _oil and C _oil are stiffness matrix and damping matrix of oil-lubricated tilted pad bearing respectively. The structural parameters of the bearing used for analysis are listed in Table 1. Similarly, the K _hub and C _hub are the force coefficient matrices of the bearing pedestal respectively. The M _hub is the mass matrix of the bearing pedestal. The whirl frequency of the shaft is expressed as a complex parameter:$$A=i\omega$$(2)

The diameter of the shaft is defined as D = 30 mm, then the support span of the rotor can be written as$$L=b\cdot D$$(3)

As shown in Figure 1, L is the parameter to be optimized. S1, S2 and S3 are the axial dimensions of the impeller, cooling device and seal, etc., respectively.

In the finite element model and the transfer matrix model, the impeller is equivalent to a rigid disk with the same mass and moment of inertia. The mass, diameter and thickness of the disc are d _i  = 100 mm, h = 10 mm and m _i  = 0.71 kg, respectively. In addition, the mass eccentricities of the disc are assumed to be e_x  = 0.5 μm and e_y  = 0.5 μm.

Besides of the thermal management for the dry gas seals and the bearings, the paper will also propose the use of a “cooling device” to add extra stiffness to support the rotor of the turbine. However, the premise is that the stiffness of the air film does not negatively affect the unbalanced response of the rotor. By “Cooling device”, it means a high-clearance journal bearings with cool CO₂ flowing through the annulus passage to cool the shaft. The cooler is not a true journal bearing. Its clearance is selected to provide adequate cooling. The flow in the annular passage in the cooling device will be simulated using computational fluid dynamics (CFD). The cooling device stiffness and damping values (as predicted in CFD) are combined with the shaft model of SCO₂ turbine. The results will give the predicted critical speed and damping for the SCO₂ turbine.

2.2 Rotor system modelling

The optimization of key parameters allows one to obtain a primitive conceptual design starting with minimal information of the object topology. Objective functions are the system's responses to these parameters, and used for defining optimization criteria. For the present work, the support system design offering the best compromise between three objective functions was chosen to be investigated on the isolated, rotating turbine rotor system. The proposed objective functions are designed to reduce the harmonic response amplitude on the rotor cantilever and to increase the distance between the rotor resonance frequency and the operating frequency. Generally, the multi-objective optimization problem can be described mathematically as$$\begin{array}{ll}\text{Minimize}\hfill & F\left(X\right)=\left(f_1\left(X\right),f_2\left(X\right),\cdots ,f_m\left(X\right)\right)\hfill \\ \text{Subject to}\hfill & X\in X^n\hfill \end{array}$$where F(X) is the m-dimensional objective function vector, f _m  (X) is the objective function of the mth dimension, X is the n-dimensional design variable, X ⁿ is an n-dimensional design space, constrained by lower and upper bounds. The parameters to be optimized along with their ranges have been summarized as$$s.t\{\begin{array}{l}3\le b\le 7\\ 0<m_{\mathrm{hub}}<m_s\\ s_l\cdot k_{oil}^w\le k_{\mathrm{hub}}\le s_u\cdot k_{oil}^w\\ s_l\cdot c_{oil}^w\le c_{\mathrm{hub}}\le s_u\cdot c_{oil}^w\\ 0<k_g<k_{oil}^w\end{array}$$(4)where $k_{oil}^w$ = 7 × 10⁸ N/m and $c_{oil}^w$ = 2 × 10⁵N . s/m are the stiffness and damping coefficient of the four-pad tilting pad bearing at working speed (50 krpm), respectively. The s _l  = 0.014 and s _u  = 10 are the scaling factor. Reserve effective space for other components while ensuring the stiffness of the shaft, so the range of parameter b is specified as [3,7]. The m_s is the mass of the shaft. The k_g is the stiffness of the gas film, assumed to be unknown.

To perform the bearing-rotor system optimization, the phase one is to establish a mathematical model of the system. To this end, two classic methods, the transfer matrix method (RTM) and the finite element method (FEM), are considered. Notably, the transfer matrix technique require fewer elements to provide an accurate system model [26]. Therefore, the RTM is used in the optimization design and the FEM is used in the simulation of rotor systems. In a RTM, the rotor is divided into several disc-shaft units, and the transfer equation between arbitrary two adjacent cross sections can be expressed as$$\left\{\begin{array}{c}\hfill f_{i+1}\hfill \\ \hfill e_{i\text{+}1}\hfill \end{array}\right\}=\left[\begin{array}{cc}\hfill u_{11}\hfill & \hfill u_{12}\hfill \\ \hfill u_{21}\hfill & \hfill u_{22}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill f_i\hfill \\ \hfill e_i\hfill \end{array}\right\}$$(5)

where, the state vectors e = ${\left[\begin{array}{cccc}\hfill X\hfill & \hfill {\theta }_x\hfill & \hfill Y\hfill & \hfill {\theta }_y\hfill \end{array}\right]}^T$ and f = ${\left[\begin{array}{cccc}\hfill M_x\hfill & \hfill Q_x\hfill & \hfill M_y\hfill & \hfill Q_y\hfill \end{array}\right]}^T$ represent the internal forces (moment and shear force) and plural trajectory of the section, respectively. The X and Y are the radial displacements of the section. The θ is the deflection angle of the section. The M and Q are moment and shear force, respectively. The u is the partitioned form of the transfer matrix corresponding to the disc-shaft units. The unbalanced response amplitude A (ω) and modal frequency ω _c are introduced by the state vector. Refer to reference [26,27] for a detailed introduction of transfer matrix.

The first objective function is defined by the maximum resonance amplitude of the cantilever disc. It is expressed as$$f_1\left(x_i\right)=A\left({\omega }_c\right)$$(6)where A = $\sqrt{X^2+Y^2}$ . The second objective function is defined by the partial derivatives of the unbalanced response function with respect to frequency, which is expected to represent the flatness of the operating frequency region if the sweep interval is sufficiently small$$f_2\left(x_i\right)=abs\left({\frac{\partial A}{\partial n}|}_{{\omega }_{work}}\right).$$(7)

In the actual numerical calculation, the second one can be written as a difference form:$$f_2\left(x_i\right)=abs\left(\frac{A_{n+1}-A_{n-1}}{2\mathrm{\backslash rm\; \Delta }\omega }\right)$$(8)where A _n+1 and A _n−1 respectively represent the amplitude of the harmonic response at the rotational speed adjacent to the working frequency, which can also be obtained from the state vector (e _i ) of the disc. Applying this transformation to the objective function can increase the discrepancy between the operating speed point and the critical speeds. The absolute value of the difference between the rotor working frequency and the nearest critical frequency is$$\mathrm{\backslash rm\; \Delta }n=abs\left({\omega }_{nc}-{\omega }_{work}\right)$$(9)

Then the third objective function is defined as$$f_3\left(x_i\right)=1/\mathrm{\backslash rm\; \Delta }n$$(10)

The formulation is employed for the evaluation of the global stiffness of a rotor system. As the stiffness of the system increases, the critical speed ratio increases, and vice versa. In some ranges of support stiffness, the variation of modal frequency is also drastic, while in other ranges it is relatively gentle [28]. These factors must be considered during optimization, so that the modal frequency of each order can be actively adjusted to a reasonable position by adjusting the rotor body or supporting parameters. The purpose of this effort is also to include further modal frequency information in objective functions. In the current problem, the optimal solutions would be the ones while all three objective functions reach their minimum.

2.3 Multi-objective optimization method

According to the above, many parameters may affect the harmonic response performance of a SCO₂ turbine rotor, including (1) hub stiffness k _hub , (2) hub damping c _hub , (3) hub mass m _hub , (4) support span b, and (5) Gas film stiffness k _g . They make up the particles in the design space, denoted as$$x_i\left(t\right)=\left[\begin{array}{cccccc}\hfill b^i\left(t\right)\hfill & \hfill m_{\mathrm{hub}}^i\left(t\right)\hfill & \hfill k_{\mathrm{hub}}^i\left(t\right)\hfill & \hfill c_{\mathrm{hub}}^i\left(t\right)\hfill & \hfill k_{gx}^i\hfill & \hfill k_{gy}^i\hfill \end{array}\right]$$(11)

Wang Peng et al. proposed the basic framework of multi-scale quantum harmonic oscillator algorithm for the first time [25]. With only one human input parameter (scale), multi-scale convergence and no training are the distinctive features of this intelligent algorithm. In the case of single-objective optimization, the update of all scales is accomplished via solving the standard deviation of the particle swarm. In this article, to treat the multi-objective functions, a non-dominated sorting algorithm (INSA) is employed to improve MQHOA. In addition, owing to the diversity of solutions for multi-objective optimization, the energy level scale matrix κ _k is defined as the standard deviation of the objective function value of the particles. The procedure of the multi-objective and −scale quantum harmonic oscillator algorithm (mMQHOA) is described below.

• Initialization

  The initial stage of the optimization procedure is to generate a set of random particles x_i and specify the convergence accuracy ( A_c ) of standard deviation. The initial scale is determined by the design space as$${\sigma }_s=x^{\max }-x^{\min }$$(12)

  Specifically, this signifies that the initial scales needs to be determined via the maximum and minimum of particles in design space.

• Energy level stabilization

  Before the start of each iteration, calculate the standard deviations of the current particles and the standard deviations of their objective function values, denoted as σ _k and κ _k . After the operation is accomplished, the routine proceeds to generate a set of new sampled particle swarm X ^N for all k particles x _i . The evolution of all particles is implemented randomly in a Gaussian distribution ($x_i^N\sim N\left(x_i,{\sigma }_s^2\right)$), which corresponds to the wave function in quantum mechanics, and serves as a constraint for solving optimization problems. The objective function values of the new particles and the current one are subsequently non-dominantly sorted by running INSA program. The optimal particle from that sort is then used to update the current particle. If the subsequent evaluation particles x ^N is superior to previously evaluated particles x _i , is replaced by x ^N . This represents an evolution of current particle swarms that are randomly distributed. The evolution is necessary as for most multi-objective and multi-dimensional engineering practices, the initial uniform sampling will not be detailed enough to identify a local extrema combination of objective functions.

  Adhering to the suggestion of Peng and Yan [25], a sub-state named energy level stabilization procedure is applied before the energy level (Energy of quantum harmonic oscillator, QHO) and scale σ _s are lowered. The sub-state is defined as the absolute value Δσ _k of the difference between the pre-evolutionary standard deviations σ _k and the post-evolutionary standard deviations ${\mathbf{\sigma }}_k^N$, i.e. $\mathrm{\Delta }{\mathbf{\sigma }}_k=abs\left({\mathbf{\sigma }}_k-{\mathbf{\sigma }}_k^N\right)$. The corresponding post-evolutionary standard deviations representing the stability of the energy level is ${\kappa }_k^N$.

  Since the evolution of each particle is independent, it is impossible to determine whether the optimal value is local or global, and the optimization will run until the energy level is stable (Δσ _k  < σ _s or ${\mathbf{\kappa }}_k^N<{\mathbf{\kappa }}_k$). The proposed strategy seeks the stable state of the objective function while balancing sampling diversity and local exploration.

• Quantum harmonic oscillator (QHO) Convergence

  After each evolution, the rank of each particle is specified in the form of non-dominated sorting. That is, fronts are ranked according to the domination concept of objective function. In order to further facilitate the energy level reduction of the design hypercube, the objective function values of the worst level are normalized within their range. It is expressed as$$f_m^{norm}\left(x_j^{wl}\right)=\frac{f\left(x_j^{wl}\right)}{\max \left(f^{wl}\right)-\min \left(f^{wl}\right)}$$(13)

  The sum of the normalized objective function for each particle is$$f_{sum}^{norm}\left(x_i^{wl}\right)={\displaystyle \sum _{m=1}^{M_o}}{f}_m^{norm}\left(x_j^{wl}\right)$$(14)The three-dimensional optimization problem is considered, M_o is set to 3. The particle with the largest sum of normalized values is the worst-evaluated particle x ^worst . The x ^worst is then replaced by the average value ${\displaystyle \overline{\mathbf{x}}}$ of all current particles, which reduces the energy level of the system while also facilitating to maintain the diversity of the evolutionary solutions. Elites are preserved for the next evolution. The optimization routine will continue to perform the procedure of energy level stabilization on this new particle swarm X . The reduction in energy level until σ _k  ≤ σ _s is termed as QHO convergence.

• Multiscale (M) Convergence

  The search accuracy depends on the standard deviation of the particle swarm, which are selected as the convergence criterions. The optimization seeks to gradually approximate the search accuracy (A _c ) by halving the scale value (σ _s  = 0.5σ _s ). The QHO procedure signifies that the algorithm explores the optimal region in the design space, while the scaling down σ _s procedure aims at further narrow reasonable regions in an effort to ensure an independence of local exploration. The scale reduction until σ _s  ≤ A _c is termed as M (Multi-scale) convergence.

  The pseudo code of the mMQHOA routine is represented in Table 2. The algorithm requires a set of initial random particles uniformly distributed in the design space to start. From the initial stage, the number of evolutionary generations depends on the dimensionality of the optimization problem and the convergence accuracy of the pre-specified standard deviation, as well as the complexity of the objective function.

  Table 2

  Pseudo code of mMQHOA.

  In case of the run does not converge within the iteration limit, the run will be terminated and random particles will be automatically generated to perform the optimization calculation again.

3 Results and discussion

3.1 Numerical experiment

For purpose of acquiring estimations of the performance and repeated fidelity, the algorithm is tested on the synthetic optimization conundrum for static multi-objective. The test equations named DTLZ2 [29] are given by$$F_D\left(X\right)=\left(f_1\left(X\right),f_2\left(X\right),f_2\left(X\right)\right)$$(15) $$\min f_1\left(X\right)=\cos \left(\frac{\pi }{2}{x}_1\right)\cos \left(\frac{\pi }{2}{x}_2\right)\left(1+g\left(X\right)\right)$$(16) $$\min f_2\left(X\right)=\cos \left(\frac{\pi }{2}{x}_1\right)\sin \left(\frac{\pi }{2}{x}_2\right)\left(1+g\left(X\right)\right)$$(17) $$\min f_3\left(X\right)=\sin \left(\frac{\pi }{2}{x}_1\right)\left(1+g\left(X\right)\right)$$(18) $$g\left(X\right)={\displaystyle \sum _{i=3}^n}{\left(x_i-0.5\right)}^2$$(19) $$s.t\{\begin{array}{l}n=\mathrm{12;};i=1,2,\cdots n\\ 0\le x_i\le 1\end{array}$$(20)

The F _D ( X ) is a multimodal concave function. The reason for using this test functions is that they represent intricate spatial surfaces. The complex surfaces with concave shapes, multiple peaks, and multiple local minima pose a moderately challenge to most optimization algorithm. They may not be able to reproduce the complex characteristics for SCO₂ bearing-rotors, but are intricate to solve and should be able to give an assessment of whether a non-dominant solutions could be obtained. Define the spatial distance as:$$s_d={\displaystyle \sum _{i=1}^{M_o}}abs\left[f_i^I\left(X\right)-f_i^P\left(X\right)\right]$$(21)

The function values of the Pareto frontier ($f_i^I\left(X\right)$) are known, so that the optimization can be evaluated as convergence when the function values lowers within a certain threshold of the fronts. The mMQHOA is run through the test problems 512 times with a fixed number of 16 stochastic particles. Each solution $f_i^P\left(X\right)$ corresponds to a Pareto front with the closest spatial distance. Therefore, the ideal Pareto front corresponding to the non-inferior optimal solution ($f_i^P\left(X\right)$) of mMQHOA is determined.

The errors of the non-dominated solutions are exemplified in Figure 3 when mMQHOA is applied to the DTLZ2 with twelve parameters. The approximation of ideal Pareto fronts by the non-dominated solutions at three different objective functions after multiple random optimizations is given. Figure 2 shows the mMQHOA's solutions (b) and the ideal Pareto fronts (a). The x-, y- and z-axis are the first-, second- and third-objective function, respectively. The errors are sorted according to the increase of the objective functions and plotted in Figure 3. It can be seen from the scatters that the randomness of Pareto fronts has been considered, and the error of mMQHOA's non-dominated solution relative to Pareto front is less than 0.15%.

Fig. 2 The ideal Pareto fronts (a) and the Pareto fronts predicted by mMQHOA (b). (a) Ideal Pareto front. (b) mMQHOA Pareto front.

Fig. 3 Error distribution of mMQHOA predicted Pareto fronts. (a) Errors of objective function 1. (b) Errors of objective function 2. (c) Errors of objective function 3.

3.2 Optimization of rotor system

3.2.1 Convergence analysis

In the multipurpose optimization that may include nonlinear force coefficients in the bearing-rotor system, somewhile it is difficult to find the optimum design scheme without utilizing optimization algorithms. Moreover, for the circumstances with conflicting objectives, acquiring a superior solution for an objective function may cause others to become worse. Therefore, for most multi-objective problems, the design space needs to be constrained and algorithm must automatically make trade-offs. The optimization was carried out in this work using the mMQHOA implemented by the in-house code. The fast and computationally program developed for this investigation is generic in nature and can be used for any other multi-objective system and rotors. For all runs and algorithms, a swarm size of sixteen is used. The calculation of the stiffness and damping of the oil-lubricated tilted pad bearings is simplified by assuming that the oil film is controlled by Reynolds equation. As the equation is still nonlinear, numerical method is used to solve it. Or introduce the infinite short bearing hypothesis and the Guimbert boundary condition to derive the explicit expression of the force coefficients [26,30]. The presence of nonlinearity complicates the optimization process.

The runs are implemented for both rotor systems of angular contact ball bearings (ACBB) and oil-lubricated tilted pad bearings (OTPB). The objective functions after 20986 (ACBB) and 33507 (OTPB) iterations shows the convergence. Figure 4 shows convergence curves for a selection of random particles. The represent the trajectory of points on the Gaussian probability of the feasible space obtained by phasing out the worst particles. In the optimization of the OTPB-rotor system, the convergences process of the objective functions are apparently more tortuous. In fact, the significant influence of the stiffness and damping of the compliant hubs is reflected in the critical speed and the amplitude of the unbalanced response, and the variation of the rotor speed will directly affect the stiffness and damping of the oil film bearing. As shown in Figure 4, this interactive and gradual running-in process is reflected in the tortuous convergence curve of objective function 3 and objective function 2.

Additionally, the performance of the mMQHOA is compared with a CIMPSOA [31], another sociology information interaction based multi-objective evolutionary algorithm. As described in Appendix A, CIMPSOA differs as it employs chaotic maps to create and update the non-linear elements in the particles. Figure 5a presents the convergence trajectory of a random routine for a OTPB-rotor system. Only the convergence trails of the global optimal particle is shown. Overall, the data shows that CIMPSO converges significantly faster than mMQHOA in multi-objective optimization. But the analysis of convergence process reveals that it seems impossible to have the minimum of all objective functions, simultaneously. For objective function 2, there is still the possibility that more iterations are needed to converge, but it is clear that the run has matured prematurely. Examining the data of the harmonic response amplitude value in Figure 5b demonstrates a wide spread of between the 2.0 and 5.5 for the dimensionless resonance amplitude. Premature convergence and especially poor repetitive fidelity imply that CIMPSOA has flaws in fine search.

Fig. 4 Comparison of convergence process of mMQHOA optimizations for the two types of bearings. (a) Objective function 1 (ACBB). (b) Objective function 1 (OTPB). (c) Objective function 2 (ACBB). (d) Objective function 2 (OTPB). (e) Objective function 3 (ACBB). (f) Objective function 3 (OTPB).

Fig. 5 Convergence procedure (a) and results (b) of CIMPSOA for optimization of OTPB-rotor system [31]. (a) Convergence procedure of CIMPSOA. (b) Harmonic response results for non-dominated solutions.

Fig. 5 Continued.

3.2.2 Harmonic response

Harmonic response analysis is performed on a rotor-bearing system to identify inherent dynamic characteristics and to generate amplitude information of the journal on critical speed and working speed. The data of harmonic response amplitude are plotted against the angular velocity of the rotor and were normalized by the eccentricity (e_x  = e_y  = 0.5 μm) of the disc on the cantilever. Figures 6 and 7 yield comparisons of the harmonic response of the disc between the empirical solution and the optimized solution for two type of bearing-rotors. The empirical particle is defined as$$X^e=\left[\begin{array}{cccccc}\hfill 4\hfill & \hfill 0.5m_s\hfill & \hfill 0.5k_{oil}^{work}\hfill & \hfill c_{oil}^{work}\hfill & \hfill \left(1e-5\right)k_{oil}^{work}\hfill & \hfill \left(1e-5\right)k_{oil}^{work}\hfill \end{array}\right]$$(22)The purpose of this design is to reduce the stiffness and reserve the damping while avoiding the influence of the gas film force on the rotor system. In the two runs shown in Figure 6, the stiffness of the hub was constrained to range A($0.1s_l\cdot k_{oil}^w\le k_{\mathrm{hub}}\le s_u\cdot k_{oil}^w$)  and range B ($s_l\cdot k_{oil}^w\le k_{\mathrm{hub}}\le s_u\cdot k_{oil}^w$) successively, so two corresponding sets of non-dominated solutions were generated respectively, as shown in Figure 6. Compared with empirical particles, there clearly is a critical frequency reduction for the optimized design, from approximately 34.33% in the hub stiffness range B to 67.16% in the range A. Figure 7 represents the harmonic response versus the angular velocity of the ACBB-rotor. Similar to the results observed in Figure 6, the amplitude is effectively attenuated at the non-dominated solutions, and the resonance amplitudes are concentrated.

Furthermore, the rotors supported by the two types of bearings have been tested by finite element method on the ANSYS software at the same design parameters conditions as imposed as parameters in the transfer matrix computations. Figure 9 shows FE results of harmonic response. The data yield that the trend witnessed in the transfer matrix models is the same as that observed in the finite element model. With the working angular velocity of 5235.99 rad/s as the constraint condition, the data shows that the maximum harmonic response amplitude of the cantilever is reduced by 71.91%, and the distance between the resonance frequency point and the working frequency point is increased by approximately 37.83%. Taking the FE results as a reference, the data obtained by the numerical simulation of FEM exhibits that the resonance frequency and resonance amplitude are reduced by approximately 4.51% (${\omega }_e^c$ for ACBB), 5.14% (${\omega }_e^c$ for OTPB) and 7.48% ($A_e^{\max }$ for ACBB), 9.18% ($A_e^{\max }$ for OTPB), respectively, and hence lower than the results predicted by transfer matrix. This shows a deviation from the smooth shape and lumped parameter models used in the numerical campaign. Nonetheless, the finite element data provides verification in both the resonance frequency and harmonic response amplitude, and can reproduce the trend with reasonable accuracy without resorting to more detailed mathematical models in the optimization algorithm.

At the current accuracy, non-inferior optimal solutions on the Pareto fronts are non-dominated capable of satisfying the requirements of rotor dynamics. It is worth noting that b, m_hub , k_hub , and c_hub reduce by 24.91%, 16.00%, 90.86%, and 92.71%, respectively, compared with empirical particles, while k_gx and k_gy enhances by approximately 3810.71% and 2641.82%, respectively, which indicates that the optimal points depends strongly on appropriate cantilever stiffness and compliant support. Several tests have noted that all the optimal points occur at b of minimum (about 3). For the practical application in which b is more necessary than other parameters. In fact, the gear and seal between the two bearings determine the lower limit of b.

Fig. 6 Harmonic response for initial design and mMQHOA optimum design of a OTPB-rotor system.

Fig. 7 Harmonic response for initial design and mMQHOA optimum design of a ACBB-rotor system.

Fig. 8 3D model used for finite element analysis.

Fig. 9 Analysis of mMQHOA Pareto fronts with finite element method.

3.2.3 Algorithm improvement

Even for OTPB with nonlinear force coefficients, mMQHOA shows stable fine search performance. Instead of a bionic interactive exploration pattern, mMQHOA uniformly initializes particles from the design space and perturbs them with multi-scale standard deviations. Further reduction of the scale results in further exploitation of the territory in the close vicinity of the current optimal solution as well as subtly probe of behavior further away from it. For each of the particles, non-dominated sorting is performed according to the objective function value predicted by the mathematical model, and the candidate particle with the highest rank is chosen as the Pareto front. In CIMPSOA, the current non-inferior optimal solutions dominate the update iteration of particles. Pertaining to the comparison of CIMPSOA and mMQHOA, the general trend across all trials indicated the interaction of information between particles resulting in accelerated convergence. And information interaction is the core of CIMPSOA. CIMPSOA performed well in convergence speed, but exhibited comparably poor results in repeated fidelity trials as well as the convergence accuracy. In Appendix A, the parameters that CIMPSOA needs to customize include inertia coefficient, self-adjustment and social adjustment weights, and the expansion factor of the chaotic map. However, mMQHOA only needs to specify the convergence accuracy, which is the advantage of updating particles with Gaussian distribution as the wave function under multi-scale. Generally speaking, the multi-scale Gaussian perturbation method performs marginally better than the search strategy of CIMPSOA, in particular for the fine optimization stage of multi- parameter and modal functions.

In sum, the evaluation of harmonic response and convergence implies that the combination of information interaction between particles and multi-scale quantum harmonic oscillator optimization algorithms can yield significant performance improvements, making the algorithm particularly attractive for fine search of Pareto fronts, such as optimal design required for rotors with nonlinear oil film forces. The following formula of information interaction inspired by CIMPSOA is used to improve the particle update iteration in the step of QHO convergence:$$x_i^{t+1}=x_i^t+0.5\left(x_{\text{opt}}^t-x_i^t\right)\text{+}0.5\left(x_{\text{mean}}^t-x_i^t\right)$$(23)where $x_{opt}^t$ and $x_{mean}^t$ are the global non-dominated particle and the mean value of all particles, in the current iteration step, respectively. The improved algorithm enables each particle to exchange information with the current non-dominated solution, in order to alternate between local fine development (Multiscale perturbation) and exploration.

Continuing from the above, in order to understand how the improved mMQHOA that introduces information interaction fare against the primary one and CIMPSOA, two independent trials are run on the rotor with OTPB and ACBB. The convergence rates for the two runs are presented in Figure 10 showing the detected three objective function values for all sixteen particles. Comparing the improved mMQHOA with CIMPSOA (as shown in Figs. 5 and 10), proves to achieve convergence with the same order of magnitude of iterative steps. For the cantilever rotor with ACBBs, the improved mMQHOA again exhibits a steady convergence fidelity, and the particle standard deviation quickly converges to a given accuracy after less than 3000 iterations. On the other hand, although the dimensionless resonance amplitude ($A_i^{\max }\left(\omega \right)$) of the information-interaction version of mMQHOA has a slightly wider distribution range (between 2.25 and 3.75) compared to the original mMQHOA (Fig. 11), the runs for OTPB-shaft systems exhibit stable and faster convergence and being able to detect the given convergence accuracy after approximately 1829 iterations. The search for the Pareto fronts are still more tortuous as a consequence of nonlinear oil film forces. In addition, compared with CIMPSOA routines that also use information interaction strategies, the improved mMQHOA performs better on repeatability fidelity as well as the convergence precision.

Fig. 10 Convergence process of OTPB rotor system optimized by improved mMQHOA. (a) Objective function 1 (ACBB). (b) Objective function 1 (OTPB). (c) Objective function 2 (ACBB). (d) Objective function 2 (OTPB). (e) Objective function 3 (ACBB). (f) Objective function 3 (OTPB).

Fig. 11 Results of improved mMQHOA for optimization of rotor system.

3.3 Sensitivity investigation

For a non-dominated solution of OTPB-rotor system shown in Table 3, parameters b and k_g are forcibly changed in order to evaluate the sensitivity of objective functions to the variation of design parameters. Hence, the variations of the resonance amplitude and its position can be observed. As shown in Figure 12, resonance amplitude always increases with increasing parameter b, but the critical frequency of the cantilever does not change significantly. Therefore, from the perspective of increasing the rigidity of the shaft, if the purpose is to achieve the lowest harmonic response amplitude, the parameter b should be set to its lowest value. However the sensitivity of the system response to gas film stiffness (k_g ) is slightly more complicated. As shown in Figure 13, the gas film stiffness is directly and positively correlated with $A_i^{\max }\left(\omega \right)$ in the median range (1E6 < k_g < 5E7), whereas in its range of high value (5E7 < k_g ), this trend is reversed. Moreover, when k_g is at a low level (k_g < 1E6), it has little effect on the unbalanced response. The sensitivity of the system to different ranges of parameters is individual, which also implies that the contradiction between the objective functions may be weakened in a certain probability.

The step changes of objective function 3 between the iterations shown in Figures 4 and 10 implies that it may have a more sensitive response to changes in design parameters than other objective functions. Figure 14 compares the effect of the presence or absence of objective function 3 on the optimization results of the rotor systems. In both cases, only the results of the respective global optimal solutions are used for comparison. As can be seen, the objective function 3 has little effect on the maximum amplitude of the harmonic response, but it directly affects the damping mode frequency. The non-dominated solution after removing the objective function 3 in the optimization routine introduces a second formant in the investigated frequency range. And the formant is near the operating speed point. Table 4 show the impact of objective function 3 on operational variables. It can be observed that the objective function 3 has the smallest impact on b, and the highest impact on the stiffness of hubs, and the damping is the second. The value of k_g is less than 1e6 N/m, so it is not considered. The stiffness of the bearing hub and cooler is the main driver for the critical speed movement, and the optimization result is that both the first and second critical frequencies avoid the operating frequency.

Fig. 12 The effects of b on the Pareto fronts.

Fig. 13 The effect of gas film stiffness on Pareto fronts.

Fig. 14 Response to unbalance for optimal solution with and without objective function 3.

4 Conclusion

This study focuses on the vibration control analysis and optimized design of a small scale SCO₂ turbine system with gas film cooler. The simulated rotor comprises a disc and a cooler, on which two bearings and hubs are designed to adjust the critical speed and resonance amplitude value. This type of bearing-rotor is used for turbine of high-inlet temperatures in the concentrated solar power industry. The approach applied is the combination of the MQHOA and the simplified fast non-dominated sorting algorithm (INSA) and the energy level stabilization method of the objective functions. Moreover, an information-interactive multi-objective MQHOA has been proposed for acceleration. The optimization then runs through an improved probe patterns: variable-scale perturbation development (local fine search) and global exploration (information interaction).

For the optimization of the SCO₂ bearing-rotor, a single operating angular velocity point was selected, and the stiffness, damping and mass of the bearing hub, the gas film stiffness of the cooler and the bearing span parameterized using six different parameters. The diameter of the rotor, the cantilever and its disk remain constant. The primary purpose of the optimization campaign is to attenuate the harmonic response amplitude while ensuring that the operating speed is not affected by the critical frequency. The conclusions are summarized as follows:

• Three different optimization runs were randomly selected and compared, where one employed the mMQHOA, one used CIMPSOA, and another improved mMQHOA. CIMPSOA converges fast but has poor repeat fidelity. Both the basic and improved mMQHOA can converge to a non-dominated solution. mMQHOA requires an average of about 27,500 design evaluations, while the improved mMQHOA requires an average of approximately 2500. Both optimization designs indicate that the harmonic response amplitude is effectively attenuated, and the resonance frequency is far away from the operating frequency. Overall, convergence rate appear to vary substantially as a consequence of the information interaction between particles, and the objective functions are predominantly decreasing in the entire convergence phases, implying that information interaction is the dominant driver for accelerating convergence.

• The harmonic response calculations of three-dimensional finite element models constructed for accurate results are performed using commercial software ANSYS. Comparing empirical and optimized design of the finite element models indicates an increase in the distance between the operating frequency and the resonance frequency of approximately 37.83% and decrease in the resonance amplitude of 71.91% confirming the trends witnessed in the observed in the transfer matrix models, and also implying that the construction of the objective function set is reasonable.

• The results show that the increase in the span of the two bearings leads to an increase in the magnitude of the unbalanced response of the disc on the cantilever. In addition, the evolution of harmonic response of the disc exhibits that the gas film stiffness of the cooler has no significant effect on the optimization results during the phase of less than 1E6 N/m. This accounts for a high stiffness produced by the brace and rotor. At median levels of air film stiffness, increasing its value attenuates the harmonic response amplitude of the cantilever, but this trend is reversed in high levels of air film stiffness.

• The analyses of harmonic response and the convergence of objective function 3 reveal that the relative position between the critical frequency and the operating frequency is most sensitive to the modification of particle elements. The stiffness of the bearing hub is the main driver for the critical speed movement. Maximizing the distance between the critical frequency and the operating frequency and minimizing the resonant amplitude may be conflicting objectives, depending on the range of design parameters, so the solutions of optimization are figured out as a set of non-dominated fronts with decision-making preferences.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work is supported by the Special Project of Industrial Cluster in National Innovation Demonstration Zone (No. 201200210400) and the Program of Henan Center for Outstanding Overseas Scientists (No. GZS2018005). The authors would like to acknowledge the support from Henan Key Laboratory for Machinery Design and Transmission System. The authors wish to acknowledge also the turbine structure and guidance provided by the team of the University of Queensland, Australia.

Appendix A: Introduction of CIMPSO

The main feature of CIMPSO (chaotic interval multi-objective particle swarm optimization algorithm) is that chaotic mapping is introduced to deal with the problem of multi-objective optimization with nonlinear elements. Based on the optimization process of the rotor system, the pseudo code of CIMPSO is shown in Table A1. The detailed procedure of this multi-objective optimization methodology can be referred to [34].

References

All Tables

All Figures

	Fig. 1 A schematic view of the investigated bearing-rotor system (a) and its cooler (b). (a) Schematic diagram of bearing-rotor system. (b) Schematic diagramc of the cooler.
In the text

	Fig. 2 The ideal Pareto fronts (a) and the Pareto fronts predicted by mMQHOA (b). (a) Ideal Pareto front. (b) mMQHOA Pareto front.
In the text

	Fig. 3 Error distribution of mMQHOA predicted Pareto fronts. (a) Errors of objective function 1. (b) Errors of objective function 2. (c) Errors of objective function 3.
In the text

	Fig. 4 Comparison of convergence process of mMQHOA optimizations for the two types of bearings. (a) Objective function 1 (ACBB). (b) Objective function 1 (OTPB). (c) Objective function 2 (ACBB). (d) Objective function 2 (OTPB). (e) Objective function 3 (ACBB). (f) Objective function 3 (OTPB).
In the text

	Fig. 5 Convergence procedure (a) and results (b) of CIMPSOA for optimization of OTPB-rotor system [31]. (a) Convergence procedure of CIMPSOA. (b) Harmonic response results for non-dominated solutions.
In the text

	Fig. 5 Continued.
In the text

	Fig. 6 Harmonic response for initial design and mMQHOA optimum design of a OTPB-rotor system.
In the text

	Fig. 7 Harmonic response for initial design and mMQHOA optimum design of a ACBB-rotor system.
In the text

	Fig. 8 3D model used for finite element analysis.
In the text

	Fig. 9 Analysis of mMQHOA Pareto fronts with finite element method.
In the text

	Fig. 10 Convergence process of OTPB rotor system optimized by improved mMQHOA. (a) Objective function 1 (ACBB). (b) Objective function 1 (OTPB). (c) Objective function 2 (ACBB). (d) Objective function 2 (OTPB). (e) Objective function 3 (ACBB). (f) Objective function 3 (OTPB).
In the text

	Fig. 11 Results of improved mMQHOA for optimization of rotor system.
In the text

	Fig. 12 The effects of b on the Pareto fronts.
In the text

	Fig. 13 The effect of gas film stiffness on Pareto fronts.
In the text

	Fig. 14 Response to unbalance for optimal solution with and without objective function 3.
In the text