Issue 
Mechanics & Industry
Volume 23, 2022



Article Number  21  
Number of page(s)  17  
DOI  https://doi.org/10.1051/meca/2022018  
Published online  02 August 2022 
Regular Article
Multiobjective optimization of a small scale SCO_{2} turbine rotor system with a shaft cooler
^{1}
School of Mechatronics Engineering, Henan University of Science and Technology, Henan 471003, China
^{2}
Henan Key Laboratory for Machinery Design and Transmission System, Henan University of Science and Technology, Luoyang, Henan 471003, China
^{3}
School of Mechanical and Mining Engineering, University of Queensland, Brisbane, Australia
^{*} email: li_jishun@163.com
Received:
24
November
2021
Accepted:
9
June
2022
The SCO_{2} turbine machines exchange energy through supercritical carbon dioxide. Their impeller has the features of hightemperature and −speed to enhance energy conversion efficiency, but the rotor needs to be cooled to be compatible with bearings and seals. The paper introduces a pivotal parameter optimization of a concentrating solar SCO_{2} turbine rotor and seeks to control the harmonic response amplitude while preserving the distance between the critical speed and the working speed. The optimization considers several parameters including bearing span, stiffness, effective mass and damping of the bearing hub, and gas film stiffness coefficients of the cooler. The optimization is accomplished using a multiobjective and −scale quantum harmonic oscillator algorithm (mMQHOA) that couples an information interaction algorithm and transfer matrix model. The application of information interaction accelerates the convergence speed of the objective functions. The verification results from the threedimensional finite element (3DFE) indicate that the nondominant design reduces resonance amplitude of the disc by approximately 71.91%, while the critical frequency increases by about 34.33% in the direction away from the operating frequency, and imply a tradeoff relationship between harmonic response amplitude and critical speed. It is further reveal that the increased gas film stiffness of cooler in the primary level interval (<1E6 N/m) has no significant effect on the harmonic response of the system. The optimization is based not only on the analysis of design parameters, but also focuses on the sensitivity of objective functions that can significantly affect dynamic performance. The models with a single variable of bearing span and film stiffness are investigated respectively, and then the sensitivity of the system response is analyzed. In addition, three different objective functions are proposed, with the purpose of constructing a universally applicable model that can be further used to optimize the analogous bearing rotor system.
Key words: Rotor system / multiobjective optimization / vibration control / Pareto front
© J. Li et al., published by EDP Sciences 2022
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Conventional steam Rankine cycle operate at a temperature of less than 550 °C [1–3], based on a tradeoff 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_{2}) energy conversion system [7]. As a potential alternative power conversion system, The SCO_{2} 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_{2} 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_{2} turbines are often extreme, a strict optimization design of the shaft bearing system is required. The integrated optimization for multiple physics including thermomechanics, rotorbearing dynamics.
Researches on these optimization algorithms for bearingrotor 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 bearingrotor 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 timeconsuming 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 multiobjective functions due to the deficiencies of the H _{∞} and H _{2} 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 multiobjective optimization algorithm in the electric machinery model to make a tradeoff 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 radialinflow 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 microwave rotor turbine efficiency. In addition to wave rotors, a hybrid of genetic algorithm with the NelderMead 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 multiobjective 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 nonparameterize 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 tradeoff 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 attentiongrabbing 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 highspeed 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 lowpower linestart synchronous motor [13]. The results show that up to 1.53% efficiency increase was obtained, in comparison to the induction motor.
In the predesign 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_{2} 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_{2} rotorbearing systems development. The present work is part of the SCO_{2} turbine system design optimization project, which is well aligned with the required innovation for highperformance turbine. The project is aimed on the rotor dynamics design optimization and vibration control of a small scale SCO_{2} turbine rotor system. Most of these previous investigations utilized singleobjective and natureinspired algorithms to solve the parameters and topology of rotorbearing system. Natureinspired algorithms (nondeterministic 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 multiobjective and multiscale 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 bearingrotor system of a small scale SCO_{2} turbine, which is shown in Figure 1, was rotordynamically 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 Oillubricated tilted pad bearings (OTPB)). Bearings and cantilever impeller are a major source of vibration in SCO_{2} 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 fourpad 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 bearingrotor system (a) and its cooler (b). (a) Schematic diagram of bearingrotor system. (b) Schematic diagramc of the cooler. 
2 Methodology
2.1 Structure of the rotorbearing 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(1)
This is the dynamic model of a series springmassdamping system. The K _{ oil } and C _{ oil } are stiffness matrix and damping matrix of oillubricated 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:(2)
The diameter of the shaft is defined as D = 30 mm, then the support span of the rotor can be written as(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 highclearance journal bearings with cool CO_{2} 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_{2} turbine. The results will give the predicted critical speed and damping for the SCO_{2} turbine.
Parameters of oillubricated tilted pad bearing.
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 multiobjective optimization problem can be described mathematically aswhere F(X) is the mdimensional objective function vector, f _{ m } (X) is the objective function of the mth dimension, X is the ndimensional design variable, X ^{ n } is an ndimensional design space, constrained by lower and upper bounds. The parameters to be optimized along with their ranges have been summarized as(4)where = 7 × 10^{8} N/m and = 2 × 10^{5}N . s/m are the stiffness and damping coefficient of the fourpad 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 bearingrotor 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 discshaft units, and the transfer equation between arbitrary two adjacent cross sections can be expressed as(5)
where, the state vectors e = and f = 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 discshaft 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(6)where A = . 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(7)
In the actual numerical calculation, the second one can be written as a difference form:(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(9)
Then the third objective function is defined as(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 Multiobjective optimization method
According to the above, many parameters may aﬀect the harmonic response performance of a SCO_{2} 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(11)
Wang Peng et al. proposed the basic framework of multiscale quantum harmonic oscillator algorithm for the first time [25]. With only one human input parameter (scale), multiscale convergence and no training are the distinctive features of this intelligent algorithm. In the case of singleobjective optimization, the update of all scales is accomplished via solving the standard deviation of the particle swarm. In this article, to treat the multiobjective functions, a nondominated sorting algorithm (INSA) is employed to improve MQHOA. In addition, owing to the diversity of solutions for multiobjective 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 multiobjective 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(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 (), 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 nondominantly 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 multiobjective and multidimensional 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 substate named energy level stabilization procedure is applied before the energy level (Energy of quantum harmonic oscillator, QHO) and scale σ _{ s } are lowered. The substate is defined as the absolute value Δσ _{ k } of the difference between the preevolutionary standard deviations σ _{ k } and the postevolutionary standard deviations , i.e. . The corresponding postevolutionary standard deviations representing the stability of the energy level is .
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 ). 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 nondominated 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(13)
The sum of the normalized objective function for each particle is(14)The threedimensional optimization problem is considered, M_{o} is set to 3. The particle with the largest sum of normalized values is the worstevaluated particle x ^{ worst }. The x ^{ worst } is then replaced by the average value 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 (Multiscale) 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 prespecified standard deviation, as well as the complexity of the objective function.
Table 2Pseudo 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 multiobjective. The test equations named DTLZ2 [29] are given by(15) (16) (17) (18) (19) (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_{2} bearingrotors, but are intricate to solve and should be able to give an assessment of whether a nondominant solutions could be obtained. Define the spatial distance as:(21)
The function values of the Pareto frontier () 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 corresponds to a Pareto front with the closest spatial distance. Therefore, the ideal Pareto front corresponding to the noninferior optimal solution () of mMQHOA is determined.
The errors of the nondominated solutions are exemplified in Figure 3 when mMQHOA is applied to the DTLZ2 with twelve parameters. The approximation of ideal Pareto fronts by the nondominated 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 zaxis are the first, second and thirdobjective 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 nondominated 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 bearingrotor 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 multiobjective problems, the design space needs to be constrained and algorithm must automatically make tradeoffs. The optimization was carried out in this work using the mMQHOA implemented by the inhouse code. The fast and computationally program developed for this investigation is generic in nature and can be used for any other multiobjective system and rotors. For all runs and algorithms, a swarm size of sixteen is used. The calculation of the stiffness and damping of the oillubricated 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 oillubricated 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 OTPBrotor 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 runningin 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 multiobjective evolutionary algorithm. As described in Appendix A, CIMPSOA differs as it employs chaotic maps to create and update the nonlinear elements in the particles. Figure 5a presents the convergence trajectory of a random routine for a OTPBrotor system. Only the convergence trails of the global optimal particle is shown. Overall, the data shows that CIMPSO converges significantly faster than mMQHOA in multiobjective 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 OTPBrotor system [31]. (a) Convergence procedure of CIMPSOA. (b) Harmonic response results for nondominated solutions. 
Fig. 5
Continued. 
3.2.2 Harmonic response
Harmonic response analysis is performed on a rotorbearing 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 bearingrotors. The empirical particle is defined as(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() and range B () successively, so two corresponding sets of nondominated 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 ACBBrotor. Similar to the results observed in Figure 6, the amplitude is effectively attenuated at the nondominated 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% ( for ACBB), 5.14% ( for OTPB) and 7.48% ( for ACBB), 9.18% ( 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, noninferior optimal solutions on the Pareto fronts are nondominated 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 OTPBrotor system. 
Fig. 7
Harmonic response for initial design and mMQHOA optimum design of a ACBBrotor 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 multiscale 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, nondominated 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 noninferior 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, selfadjustment 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 multiscale. Generally speaking, the multiscale 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 multiscale 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:(23)where and are the global nondominated 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 nondominated 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 () of the informationinteraction 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 OTPBshaft 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 nondominated solution of OTPBrotor 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 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 nondominated 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.
Comparison of Pareto fronts and empirical solutions.
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. 
The impact of objective function 3 on the Pareto fronts.
4 Conclusion
This study focuses on the vibration control analysis and optimized design of a small scale SCO_{2} 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 bearingrotor is used for turbine of highinlet temperatures in the concentrated solar power industry. The approach applied is the combination of the MQHOA and the simplified fast nondominated sorting algorithm (INSA) and the energy level stabilization method of the objective functions. Moreover, an informationinteractive multiobjective MQHOA has been proposed for acceleration. The optimization then runs through an improved probe patterns: variablescale perturbation development (local fine search) and global exploration (information interaction).
For the optimization of the SCO_{2} bearingrotor, 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 nondominated 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 threedimensional 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 nondominated fronts with decisionmaking 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 multiobjective particle swarm optimization algorithm) is that chaotic mapping is introduced to deal with the problem of multiobjective 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 multiobjective optimization methodology can be referred to [34].
References
 Y. Ahn, S.J. Bae, M. Kim, S.K. Cho, S. Baik, J.I. Lee, J.E. Cha, Review of supercritical CO_{2} power cycle technology and current status of research and development, Nucl. Eng. Technol. 47 , 647–661 (2015) [Google Scholar]
 F. Crespi, G. Gavagnin, D. Sanchez, G.S. Martinez, Supercritical carbon dioxide cycles for power generation: a review, Appl. Energy 195 , 152–183 (2017) [CrossRef] [Google Scholar]
 Y.P. Liu, Y. Wang, D.G. Huang, Supercritical CO_{2} brayton cycle: a stateoftheart review, Energy 189 (2019) [Google Scholar]
 F. Crespi, D. Sanchez, G.S. Martinez, T. SanchezLencero, F. JimenezEspadafor, Potential of supercritical carbon dioxide power cycles to reduce the levelised cost of electricity of contemporary concentrated solar power plants, Appl. Sci.Basel. 10 (2020) [Google Scholar]
 P. Dutta, High temperature solar receiver and thermal storage systems, Appl. Therm. Eng. 124 , 624–632 (2017) [CrossRef] [Google Scholar]
 J.M. Yin, Q.Y. Zheng, Z.R. Peng, X.R. Zhang, Review of supercritical CO_{2} power cycles integrated with CSP, Int. J. Energy Res. 44 , 1337–1369 (2020) [CrossRef] [Google Scholar]
 M.M. Ehsan, Z.Q. Guan, H. Gurgenci, A. Klimenko, Feasibility of dry cooling in supercritical CO_{2} power cycle in concentrated solar power application: review and a case study, Renew. Sustain. Energy Rev. 132 (2020) [Google Scholar]
 M.T. White, G. Bianchi, L. Chai, S.A. Tassou, A.I. Sayma, Review of supercritical CO_{2} technologies and systems for power generation, Appl. Thermal Eng. 185 (2021) [Google Scholar]
 R. Rzadkowski, G. Zywica, T.Z. Kaczmarczyk, A. Koprowski, K. Dominiczak, R. Szczepanik, M. Kowalski, Design and investigation of a partial admission radial 2.5kW organic Rankine cycle microturbine, Int. J. Energy Res. 44 , 11029–11043 (2020) [CrossRef] [Google Scholar]
 H. Heidari, P. Safarpour, Optimal design of support parameters for minimum force transmissibility of a flexible rotor based on H and H2 optimization methods, Eng. Optim. 50 , 671–683 (2018) [CrossRef] [MathSciNet] [Google Scholar]
 J. Lim, S. Shin, Y. Kee, Optimization of rotor structural design in compound rotorcraft with lift offset, J. Am. Helicopter Soc. 61 (2016) [Google Scholar]
 O. Laldin, S.D. Sudhoff, S. Pekarek, An analytical design model for wound rotor synchronous machines, IEEE Trans. Energy Convers. 30 , 1299–1309 (2015) [CrossRef] [Google Scholar]
 L. Knypinski, K. Paweloszek, Y. Le Menach, Optimization of lowpower linestart PM motor using gray wolf metaheuristic algorithm, Energies 13 (2020) [Google Scholar]
 A.L.E. Sarmiento, R.G.R. Camacho, W. de Oliveira, E.I.G. Velasquez, M. Murthi, N.J.D. Gautier, Design and offdesign performance improvement of a radialinflow turbine for ORC applications using metamodels and genetic algorithm optimization, Appl. Therm. Eng. 183 (2021) [Google Scholar]
 S. Tuchler, C.D. Copeland, Numerical optimisation of a microwave rotor turbine using a quasitwodimensional CFD model and a hybrid algorithm, Shock Waves 31 , 271–300 (2021) [CrossRef] [Google Scholar]
 L. Witanowski, P. Klonowicz, P. Lampart, T. Suchocki, L. Jedrzejewski, D. Zaniewski, P. Klimaszewski, Optimization of an axial turbine for a small scale ORC waste heat recovery system, Energy 205 (2020) [Google Scholar]
 K.J. Shibu, K. Shankar, C.K. Babu, G.K. Degaonkar, Multiobjective optimisation of a small aircraft turbine engine rotor system with selfupdating Rayleigh damping model and frequencydependent bearingpedestal model, Proc. Inst. Mech. Eng. C 233 , 5710–5723 (2019) [CrossRef] [Google Scholar]
 S. Hiruma, M. Ohtani, S. Soma, Y. Kubota, H. Igarashi, Novel hybridization of parameter and topology optimizations: application to permanent magnet motor, IEEE Trans. Magn. 57 (2021) [CrossRef] [Google Scholar]
 A. Nag, H. Ramachandran, A. Shriwastava, Optimization of the interference parameters of an Orbit motor using genetic algorithm, Proc. Inst. Mech. Eng. C 234 , 4478–4492 (2020) [CrossRef] [Google Scholar]
 T. Ishikawa, K. Nakayama, N. Kurita, F.P. Dawson, Optimization of rotor topology in PM synchronous motors by genetic algorithm considering cluster of materials and cleaning procedure, IEEE Trans. Magn. 50 (2014) [Google Scholar]
 W. Kwak, Y. Lee, Optimal design and experimental verification of piezoelectric energy harvester with fractal structure, Appl. Energy 282 (2021) [Google Scholar]
 J.M. Ahn, J.C. Son, D.K. Lim, Optimal design of outerrotor surface mounted permanent magnet synchronous motor for cogging torque reduction using territory particle swarm optimization, J. Electr. Eng. Technol. 16 , 429–436 (2021) [CrossRef] [Google Scholar]
 R.R. Mutra, J. Srinivas, R. Rzadkowski, An optimal parameter identification approach in foil bearing supported highspeed turbocharger rotor system, Arch. Appl. Mech. 91 , 1557–1575 (2021) [CrossRef] [Google Scholar]
 G.F. Gomes, J.A.S. Chaves, F.A. de Almeida, An inverse damage location problem applied to AS350 rotor blades using bat optimization algorithm and multiaxial vibration data, Mech. Syst. Signal Process. 145 (2020) [Google Scholar]
 W. Peng, H. Yan, L. Bo, X. Qianhe, Multiscale quantum harmonic oscillator optimization algorithm (Posts and Telecom Press, Beijing, 2016) [Google Scholar]
 W. Zheng, Rotor Dynamics Design of Rotating Machinery (Tsinghua University Press, Beijing, 2015) [Google Scholar]
 R. Albuquerque, D.L. Barbosa, Evaluation of bending critical speeds of hydrogenerator shaft lines using the transfer matrix method, Proc. Inst. Mech. Eng. C 227 , 2010–2022 (2013) [CrossRef] [Google Scholar]
 L. Ming, L. Zigang, Nonlinear Vibration of RotorBearing System under Holonomic Constraints (Science Press, Beijing, 2016) [Google Scholar]
 C.C. Coello, G.B. Lamont, D.A.v. Veldhuizen, Evolutionary Algorithms for Solving MultiObjective Problems (Springer, Boston, MA, New York, 2007) [Google Scholar]
 P. Oscar, S. Beno, Theory of Hydrodynamic Lubrication (McGrawHill Book Company, NewYork, 1961) [Google Scholar]
 J. Li, H. Gurgenci, J. Li, L. Li, Z. Guan, F. Yang, Optimal design to control rotor shaft vibrations and thermal management on a supercritical CO_{2} microturbine, Mech. Ind. 22 (2021) [Google Scholar]
Cite this article as: J. Li, H. Gurgenci, Z. Guan, J. Li, L. Li, Y. Xue, Multiobjective optimization of a small scale SCO_{2} turbine rotor system with a shaft cooler, Mechanics & Industry 23, 21 (2022)
All Tables
All Figures
Fig. 1
A schematic view of the investigated bearingrotor system (a) and its cooler (b). (a) Schematic diagram of bearingrotor 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 OTPBrotor system [31]. (a) Convergence procedure of CIMPSOA. (b) Harmonic response results for nondominated solutions. 

In the text 
Fig. 5
Continued. 

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

In the text 
Fig. 7
Harmonic response for initial design and mMQHOA optimum design of a ACBBrotor 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.