© J. Li et al., Published by EDP Sciences 2021

1 Introduction

Supercritical CO2 Brayton cycle is proposed as the next-generation technology in a diverse range of applications including nuclear power generation, fossil fuel power generation, and concentrated solar power [1,2]. Thermodynamic energy conversion system is the key issue in the future power generation. As described in the relevant literature [3], under the condition that the inlet temperature of the turbine is greater than 550  ° C, the ultra-supercritical (USC) steam cycle steam is demanded to improve the efficiency of the traditional Rankine cycle which has been widely utilized over the past decades, but the USC may interact with structural materials to create a potential threat to the circulation system.As the turbine inlet temperature gets higher, the steam cycle in any form becomes even less suitable. Compared with steam at the same temperature, the SCO2 has the advantages of less corrosion and less work required for compression. As a result, the SCO2 Brayton cycle can be considered as an alternative power conversion system, which can operate at high inlet turbine temperatures in a safe and reliable fashion [3,4].

The CO2 fluid becomes approximately incompressible near the critical point (30.98  ° C and 7.38 MPa), allowing it has a very high power density in supercritical state leading to compact footprint, small turbomachinery and micro-channel heat exchangers are capitalized in SCO2 cycles [5]. There has been significant research focus on SCO2 Brayton cycle for the next generation nuclear reactor systems and other power plants. Most of the investigations have been theoretical including thermodynamic cycle optimization directed towards determination of the optimum cycle lay-outs in different applications. Some higher fidelity numerical simulations using computational fluid dynamic (CFD) have been conducted modelling components of Brayton cycles, such as heat conduction medium, turbine, and heat receiver [6–10].

Recent interest in supercritical CO2 cycles has been mostly directed at large-scale power generation using utility-scale axial turbines, e.g. in coal fired power generation with carbon capture and nuclear power generation. SCO2 cycles can also enable modular scalable and efficient Concentrating Solar Thermal (CST) power generation with storage [11–13]. Average cost of electricity from combined CST + PV solar power generators would be cheaper than diesel generation and would offer day and night affordable electricity to communities around the world unconnected to national power grids. The focus of this paper is on SCO2 radial-in-flow turbines suitable for small-scale power generation at turbine sizes from 0.3 to 30 MWe. As mentioned above, the immediate market for this type of turbine is off-grid power generation using CST with thermal storage.

As shown in Figure 1 our proposed radial-in-flow turbine configuration for a SCO2 cycle employs an overhung shaft supported by two bearings, which use conventional technology (hydrodynamic or magnetic bearings). The reason for using an overhung shaft is to protect the two off-the-shelf bearings against the high-temperature high-pressure cycle fluid using at least two seals. This overhung configuration coupled with very high rotational speeds pose a significant challenge in terms of shaft vibrations. The University of Queensland and Henan University of Science and Technology have a joint project investigating design options to control vibrations and rotor Thermal conduction management for this configuration. The following options are being considered for the two radial bearings:

• Magnetic bearings with controllable support stiffness and damping

• Oil-lubricated tilted pad bearings in rigid hubs

• Oil-lubricated tilted pad bearings in compliant hubs

At these high rotational speeds, the expected life in operating hours will be limited for rolling-element bearings. Therefore, they are not considered in this study. The analysis will be conducted for a 300-kWe rotor supported by Oil-lubricated tilted pad bearings in compliant hubs rotating at a speed of 50000 RPM. The low end of the power range is used because small SCO2 turbines offer the greatest challenge in terms of shaft dynamics due to their high rotational speed. The optimum shaft speed gets lower as the turbines get larger [8].

For the purpose of this study, the rotor geometry is represented as shown in Figure 2. The lengths of the labyrinth seal, the cooling zone, the DGS are constrained by other requirements. For example, the length of the labyrinth seal is constrained by leakage considerations; the cooling zone needs to be long enough to cool the rotor shaft to a temperature the DGS can tolerate; the length of the DGS section is determined by the seal manufacturer's design. The values given in the table are reasonable approximations of the actual values. The bearing spacing is also expressed as bD and the effect of varying b on shaft dynamics will be examined as the main constraint on bearing spacing is placed by shaft dynamics, which is the topic of the present study. The cooling device is used to prevent blade temperatures from affecting the rotor support system and to improve cantilever stiffness as a gas bearing. The necessary first step in cooling device design is the prediction of the dynamic characteristics: based on the harmonic response of rotor dynamics, the vibration mass, stiffness and damping of the compliant hubs, and stiffness coefficient of the cooling device are optimized by multi-objective numerical calculation. The design of the cooling device was executed by employing a three-dimensional CFD-based methodology while the performance and efficiency of the cooling device were investigated.

Fig. 1 Supercritical CO2 radial-in-flow turbine lay-out.

Fig. 2 Approximate rotor shaft geometry.

2 Optimal design of rotor system

2.1 Rotor system structure

In the turbine rotor system shown in Figure 2, the composite stiffness and damping of the main support with a distance of L can be expressed as [14]$$K_e=\frac{\left[K_{oil}+A{C}_{oil}\right]\left[K_{hub}+A{C}_{hub}+A^2{m}_{hub}\right]}{\left[K_{oil}+A{C}_{oil}+K_{hub}+A{C}_{hub}+A^2{m}_{hub}\right]}$$(1)

For tilting pad bearing support system with flexible bearing pedestal, _{K oil} and C _oil are the stiffness and damping matrix of the bearing respectively. K _hub and C _hub are the stiffness and damping matrix of the bearing pedestal respectively. The whirl frequency of rotor is expressed as complex number$$A=i\omega$$(2)

The diameter of the rotor is D, and the bearing spacing is defined as$$L=b\cdot D$$(3)

where parameter b is the coefficient to be optimized. S1, S2 and S3 are axial dimensions of impeller, cooling device and sealing mechanism respectively. The lubricating medium for cooling deviceis SCO₂. As an auxiliary support, its stiffness is a nonlinear function of rotor speed (n) and rotor displacement (s), which can be expressed as$${{\displaystyle \bar{k}}}_g=k\left(n,s\right)$$(4)

In order to establish the finite element and transfer matrix model, the impeller is equivalent to a rigid disc with diameter d = 100 mm and thickness h = 10 mm, and the mass of the disc is 0.71 kg. The reason for this equivalent operation is that the disc and impeller have the same diameter and pole moment of inertia.

2.2 Objective function of rotor system optimization

The cooling and supporting parameters of the shafting of the SCO₂ solar generator directly affect its critical speed and the amplitude of the imbalanced response. Therefore, the bearing spacing L, the stiffness k_hub, the damping c_hub and the mass m_hub of the bearing pedestal, and the stiffness ${{\displaystyle \bar{k}}}_g$ of the cooling device are taken as the input variables of the system to be optimized. The mathematical model of rotor system can be established by the classical Riccati transfer matrix method [14]. The rotor is divided into several disc-axis elements, and the relation between any two adjacent cross sections can be described 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 output parameters e and f are the internal forces (bending moment and shear force) of the plural trajectory and section respectively. The u is the transfer matrix corresponding to the disc-axis element. A detailed introduction of the transfer matrix can be found in reference [14]. Rotor system design is a multi-objective optimization problem, which can be used as fitness evaluation functions including.

• The maximum peak of the unbalance response of the rotor

$$A\left(t\right)=\max \left(\begin{array}{cccc}\hfill A_{peak}^1\hfill & \hfill A_{peak}^2\hfill & \hfill \cdots \hfill & \hfill A_{peak}^i\hfill \end{array}\right)$$(6)

The $A_{peak}^i$ in the formula can be obtained from the complex parameter e _i of the current node.

• The absolute value of the difference between the operating speed and the nearest critical speed

$$\mathrm{\Delta }n\left(t\right)=\text{min}\left(abs\left(n_{\text{critical}}-n_{\text{w}\text{o}\text{r}\text{k}}\right)\right)$$(7)

Where, n_critical is the critical frequency closest to the operating frequency n_work.

• The derivative of the amplitude of the unbalance response at the working speed point to the frequency

$$A_{DC}\left(t\right)=\max \left(-abs\left(\frac{\partial A}{\partial n}|\begin{array}{c}\hfill \mathrm{}\hfill \\ \hfill n_{\text{w}\text{o}\text{r}\text{k}}\hfill \end{array}\right)\right)$$(8)

In the actual numerical calculation, it can be written as a difference form:$$A_{DC}\left(x\right)=abs\left(\frac{A_{n+1}-A_{n-1}}{2\mathrm{\Delta }n}\right)$$

where, A_n+1 and A_n−1 respectively represent the amplitude of the unbalanced response at the rotational speed point adjacent to the operating frequency point, which can also be obtained from the complex vibration parameter e _i.

The fitness function of constraint optimization can be written as$$f\left(t\right)=f\left(A\left(t\right),\mathrm{\Delta }n\left(t\right),A_{DC}\left(t\right)\right)$$(9)

A certain space is reserved for the output part while ensuring the stiffness of the rotor. Therefore, the constraint range of parameter b is$$s.t\begin{array}{c}\hfill \mathrm{}\hfill \end{array}3\le t\le 6$$(10)

The constraint range of relevant parameters of bearing pedestal can be taken as$$s.t\begin{array}{c}\hfill \mathrm{}\hfill \end{array}\{\begin{array}{c}\hfill m_{hub}>0\hfill \\ \hfill k_{hub}\ge s\cdot k_{oil}^w\hfill \\ \hfill c_{hub}\ge s\cdot c_{oil}^w\hfill \end{array}$$(11)

where, $k_{oil}^w$ and $c_{oil}^w$ are respectively the Maximum stiffness and Maximum damping of the tilting tile bearing at working speed, and s is the scaling factor. The stiffness of SCO2 cooling device is a nonlinear function of rotor speed and displacement. In practical application, the stiffness of gas bearing can be constrained as$$s.t_0<{{\displaystyle \bar{k}}}_g\le K_{upper}$$(12)

where, K_upper is the upper limit of gas bearing stiffness in practical application, which is determined by the pressure of lubricating medium and design parameters ofbearings. In rotor dynamics, for every harmonic response analysis, the stiffness of gas bearing brought into the system is a two-dimensional matrix related to speed and rotor displacement. Due to the nonlinearity of gas film stiffness, the stiffness of cooling device is unknown during numerical analysis. This brings challenges to the application of conventional genetic algorithm, simulated annealing algorithm and particle swarm optimization algorithm in the optimization of shafting parameters. In order to deal with the above problems, a chaos interval multi-objective particle swarm optimization algorithm (CIMPSO) was proposed, and the rotor system parameters of the turbine are optimized.

Within the constraint range of stiffness of the cooling device, the characteristic number k_ge is taken as an element in the particle vector, and the chaos interval particle to be optimized can be described as$$X_i\left(t\right)=\left[\begin{array}{ccccc}\hfill b^i\left(t\right)\hfill & \hfill m_{hub}^i\left(t\right)\hfill & \hfill k_{hub}^i\left(t\right)\hfill & \hfill c_{hub}^i\left(t\right)\hfill & \hfill k_{ge}^i\hfill \end{array}\right]$$(13)

2.3 Optimization of rotor system based on CIMPSO

Particle swarm optimization (PSO) is a kind of computing technology based on self-cognition and social psychology. Every potential feasible solution in optimization space is abstracted as a particle vector X _i. PSO starts with a random initial population in the constrained space, and each particle is defined by velocity and displacement. In the space, individual particles achieve information sharing mechanism while local optimization. The correct behavior of the particles is strengthened, while being imitated by other particles, and finally converges to the optimal solution or the non-inferior optimal target domain (Pareto front). The standard particle swarm optimization algorithm uses the following formula to iteratively update particles [15]:$$V_i\left(t+1\right)=w_v{V}_i\left(t\right)+c_1{r}_1\left[P_i\left(t\right)-X_i\left(t\right)\right]+c_2{r}_2\left[P_g\left(t\right)-X_i\left(t\right)\right]$$(14) $$X_i\left(t+1\right)=w_x{X}_i\left(t\right)+V_i\left(t+1\right)$$(15)

Suppose the number of initialized particles is N, i = 1, 2..., N, means the sequence number of particles. The fitness of each particle is determined by the objective function f (t). In the iterative process, each particle remembers and tracks the current individual optimal particle P _i (t) and the global optimal particle P _g (t). w_v and w_x are inertial weight and forgetting inertial weight, which determine the global optimization ability of the algorithm.r₁ and r₂ are random numbers between 0 and 1, indicating the cognitive ability of the particle.c₁ and c₂ is the learning factor, indicating the learning ability of particles.

For the shafting in this paper, the stiffness of the cooling device is a two-dimensional nonlinear elements (DNE) that varies according to the rotation speed and displacement of the shaft. The values of stiffness elements are constrained by applications. Set the upper and lower limits of the constraints as$$\begin{array}{l}\varphi \left(t\right)=\left(1-\alpha \right)k_{ge}\left(t\right)\\ \sigma \left(t\right)=\left(1+\beta \right)k_{ge}\left(t\right)\end{array}$$(16)

where, α and β are expansion coefficients based on the characteristic number of particle elements. In order to achieve particle swarm optimization with DNE, chaos variables are introduced here. The logistic equation of a typical chaotic system [15]$$z_{n+1}=\mu z_n\left(1-z_n\right),n=\mathrm{1,2},\cdots N$$(17)can be used to obtain random chaotic variables. Chaotic variables have the characteristics of regularity, ergodicity and randomness in a certain range. Where, μ is the control parameter, z_n ∈ (0, 1). After updating the particle X _i at each iteration step, the stiffness characteristic number k_ge of the cooling device is mapped to the Logistic equation to initialize the chaotic variables$$z_1\left(t\right)=\frac{k_{ge}\left(t\right)-\varphi \left(t\right)}{\sigma \left(t\right)-\varphi \left(t\right)}$$(18)where, t represents the current iteration step. Substitute z₁ (t) into Logistic equation for iteration to produce chaotic variable sequence z_i+1 (t) (i = 1, 2, ⋯ N). The chaotic sequence z_i+1 (t) is inversely mapped to the solution space to obtain$$k_i^{chaos}\left(t\right)=\varphi \left(t\right)+\left[\sigma \left(t\right)-\varphi \left(t\right)\right]z_i\left(t\right),\left(i=\mathrm{1,2},\cdots ,N\right)$$(19)

$k_i^{chaos}$ is used to replace the stiffness characteristic number in the current particle swarm, and a new input particle with chaotic characteristic is obtained$$X_i^{chaos}\left(t\right)=\left[\begin{array}{ccccc}\hfill b\left(t\right)\hfill & \hfill m_{hub}\left(t\right)\hfill & \hfill k_{hub}\left(t\right)\hfill & \hfill c_{hub}\left(t\right)\hfill & \hfill k_i^{chaos}\left(t\right)\hfill \end{array}\right]$$(20)

where, i = 1, 2, ... , N represents the number of particles. The above mentioned iterative process of generating chaotic sequences from variables is called chaotic mapping, denoted as$$Chaos\left(\text{.}\right)$$(21)

Particle swarm $X_i^{chaos}\left(t\right)$ with chaotic characteristics represents all possible stiffness coefficients determined by stiffness eigenvalues. This approach transforms the spatial complexity of the system into time complexity and then solves the non-inferior optimal target domain. The fitness of chaotic characteristic particle $X_i^{chaos}\left(t\right)$ is calculated in solution space. Each particle was evaluated simultaneously with three fitness values A (t), Δn (t) and A_DC (t), corresponding to three individual historical optimal particles $P_i^A\left(t\right)$, $P_i^{\mathrm{\Delta }n}\left(t\right)$ and $P_i^{ADC}\left(t\right)$ respectively. Then, perform fast non-dominated sorting is performed on the fitness value, and the global optimal particle is selected from the first-order particles, denoted as$$P_g\left(t\right)=\text{M}fit\left[P_i^{\mathrm{}}\left(t\right)\right]$$(22)

where, M fit (⋅) refers to the particle with the Comprehensive optimal fitness in the current iteration. Sub-objective function normalization:$$M_A^i=\frac{A^i\left(x\right)-\text{min}\left(A\right)}{\max \left(A\right)-\text{min}\left(A\right)}$$(23) $$M_{\mathrm{\Delta }n}^i=\frac{1/\mathrm{\Delta }{n}^i\left(x\right)-\text{min}\left(1/\mathrm{\Delta }n\right)}{\max \left(1/\mathrm{\Delta }n\right)-\text{min}\left(1/\mathrm{\Delta }n\right)}$$(24) $$M_{ADC}^i=\frac{A_{DC}^i\left(x\right)-\text{min}\left(A_{DC}\right)}{\max \left(A_{DC}\right)-\text{min}\left(A_{DC}\right)}$$(25)

The comprehensive fitness value of each individual particle was obtained by$$f_i\left(x\right)={\varphi }_1\cdot M_A^i+{\varphi }_2\cdot M_{\mathrm{\Delta }n}^i+{\varphi }_3\cdot M_{ADC}^i$$(26)Where, ϕ_k represent the optimal particle weight factor, and ϕ₁ + ϕ₂ + ϕ₃ = 1, i is the number of particle.

In the multi-objective optimization problem of rotor system, each fitness function may be gamed against each other, so every time the particle fitness is evaluated, the optimal solution for degradation needs to be eliminated. The particle $X_{worst}^A\left(t\right)$, $X_{worst}^{\mathrm{\Delta }n}\left(t\right)$ and $X_{worst}^{ADC}\left(t\right)$ with the worst fitness values of A (t), Δn (t) and A_DC (t) in the current iteration step are respectively searched. Three new random particles are obtained by chaos mapping the global optimal particle P _g (t), and the worst particle is replaced$$\begin{array}{l}{X}_{worst}^A\left(t\right)=Chaos\left[P_g\left(t\right)\right]\\ X_{worst}^{\mathrm{\Delta }n}\left(t\right)=Chaos\left[P_g\left(t\right)\right]\\ X_{worst}^{ADC}\left(t\right)=Chaos\left[P_g\left(t\right)\right]\end{array}$$(27)

In the current iteration step, k_ge (t) is still used as the element of the stiffness characteristic number of the replaced particles. Finally, new particles are iteratively obtained according to standard particle swarm optimization (Eqs. (14) and (15)) algorithm:$$X_i^{chaos}\left(t+1\right)=\left[\begin{array}{ccccc}\hfill b\left(t+1\right)\hfill & \hfill m_{hub}\left(t+1\right)\hfill & \hfill k_{hub}\left(t+1\right)\hfill & \hfill c_{hub}\left(t+1\right)\hfill & \hfill k_i^{chaos}\left(t+1\right)\hfill \end{array}\right]$$(28)

The optimization method described above is called chaotic intervalmulti-objective particle swarm optimization algorithm (CIMPSO). Compared with the standard chaotic particle swarm optimization (ChPSO) algorithm, the CIMPSO proposed in this paper can perform spatial optimization for multi-objective particle swarm models with two-dimensional nonlinear elements. When individual and global optimal particles do not evolve with iterations, then each particle of the particle swarm updates independently, and finally all particles tend to the global optimal solution or non-inferior optimal solution. However, all particles in the particle swarm contain chaotic maps with Characteristic number elements, which ensures effective optimization calculation of the interval of chaos interval elements within the range of chaotic maps. The convergence test for CIMPSO is in Appendix 1.

2.4 Induction of CIMPSO algorithm

The linearized oil film force model of oil-lubricated bearings is suitable for rotor micro-vibration. Some rotating machinery has the characteristics of high speed, lowflexibility and large span, however. The linear stiffness model based on small disturbance may have a large deviation from the actual situation. The rotor may exhibit nonlinear dynamic behavior under the action of nonlinear oil film force. In the process of dealing with the optimization problem of high-speed rotor systems, the main difficulty faced is the nonlinearity of the parameters such as the stiffness of the oil-lubricated bearings and the cooling device. The main advantage of CIMPSO isthat chaotic mapping is introduced to solve the problem of multi-objective optimization with nonlinear elements, the pseudo code of the algorithm is shown in Table 1. Based on the optimization process of the rotor system, the algorithm is mainly divided into the following steps.

Step 1: Initialize relevant parameters, including learning factor c₁ and c₂, particle swarm size N, optimization times T, velocity inertia weight w_v and position inertia weight w_x, expansion coefficient α and β, control parameters μ.

Step 2: All particle positions are randomly initialized within the constraint range. The position of the particle is denoted as X _i (t, j)(j represents the position number of the nonlinear element in the particle vector X, t represents the number of iteration steps). Determine the characteristic number X _i (t, c) of the nonlinear element in the particle vector. The X _i (t, c) is substituted into the Equation (21) to generate the chaotic variable sequenceand replace the characteristic number elements in all particle vectors:$$X_i\left(t,c\right)=Chaos\left[X_i\left(t,c\right)\right]$$

Step 3: Calculate the fitness of the particle vector. Each fitness corresponds to a historical individual optimal value of the particle, denoted as $P_i^1$,$P_i^2$,…$P_i^k$ (Assume that the number of objective functions is k). Perform non-dominated sorting on all individual optimal values of the current iteration step to obtain the first dominance layer particle swarm.$$P_d=\left[P_1\left(t\right),P_2\left(t\right),\cdots ,P_s\left(t\right)\right]$$

According to the operator's preference, the particle with the best fitness is selected from the first dominating layer as the global optimal particle.$$P_g=\max fit\left(P_d\right)$$

Step 4: Generate a sequence of chaotic variables:$$X_k^{chaos}\left(t\right)=Chaos\left[P_g\left(t\right)\right]$$where, k is the number of the objective function. Respectively replace the particles with the worst fitness that the decision maker considers in the current iteration step with the chaotic variable $X_k^{chaos}\left(t\right)$. The characteristic number element of each particle remains unchanged. Finally, the particle vector is updated iteratively according to the standard particle swarm optimization algorithm equations (14) and (15).

2.5 Optimal calculation

As mentioned above, the parameters to be optimized for the rotor system of the SCO₂ turbine include the bearing spacing L = Db of the main supporting system, the mass m_hub, stiffness k_hub and damping c_hub of bearing pedestal, and the stiffness ${{\displaystyle \bar{k}}}_g$ of cooling device. Where ${{\displaystyle \bar{k}}}_g$ is a chaos interval element, namely anonlinear function of rotor speed and displacement, and its characteristic number is k_ge. The constraint condition of each element is$$s.t\{\begin{array}{c}\hfill 3\le b\le 6\hfill \\ \hfill 0<m_{hub}\hfill \\ \hfill k_{hub}\ge s_1\cdot k_{oil}^w\hfill \\ \hfill c_{hub}\ge s_2\cdot c_{oil}^w\hfill \\ \hfill 0<k_{ge}<K_{upper}\hfill \end{array}$$(29)where,$k_{oil}^w=7\times {10}^8\text{N/m}$ and $c_{oil}^w=2\times {10}^5\text{N}\cdot \text{s/m}$ are the stiffness and damping coefficient of the tilting pad bearing at working speed respectively. And, s₁ = s₂ = 0.8,K_upper = 1 × 10⁷ N/m. A chaotic interval multi-objective particle swarm model was established, and parameter initialization was shown in Table 2. The particles were randomly initialized in the constraint space, and five groups of results were randomly selected after repeated optimization calculations, as shown in Table 3. The average particle ${\displaystyle \bar{x}}$ is obtained from the element average of all the optimized particles.

With reference to design experience, the bearing spacing is 3–5 times the diameter of the rotator, and the distance here is 3. The vibration mass of the bearing pedestal is half of the rotor mass (2.49 kg), In this analysis, the mass is 1.2 kg; The stiffness of the bearing hubs is half of the tilting pad bearing's, and the damping coefficient is equal to that of the tilting pad bearing; The stiffness of the gas bearing takes the middle value of the constraint space. Then the empirical particle parameter is$$x_e=\left[\begin{array}{ccccc}\hfill 3\hfill & \hfill 1.25\hfill & \hfill k_{oil}^w/2\hfill & \hfill c_{oil}^w\hfill & \hfill K_{upper}/2\hfill \end{array}\right]$$(30)

The particles were randomly initialized in the constraint space, and optimized for 25 times consecutively. Each time, the objective function was evaluated 200 times, and 25 groups of non-inferior optimal solutions were obtained. The unbalanced response generated by bringing the 25 groups of non-inferior optimal solutions into the rotor system are shown as the solid lines in the Figure 3. The dot-dash line shown in the figure represents the unbalanced response of the rotor system with empirical particle X _e, and the imbalanced response of the rotor calculated with the non-inferior optimal particles are shown as the solid lines. It can be seen from the figure that the relative peak value before optimization is 13.73 − 1.94 = 11.79. The imbalanced response of the rotor calculated with the non-inferior optimal particles in Table 3 is shown as the solid line in Figure 3. Compared with the dashed-dotted line, the solid line shows a minimum peak attenuation rate of 68.3%, and the response curve is more gentle, which reflects the robustness of CIMPSO algorithm. The difference between the working speed (n_work) and the nearest critical speed accounts for about 30% of the n_work, which ensures that the rotor works in a smooth steady whirling region.

Fig. 3 Rotor imbalance response before and after optimization design.

3 Design of cooling device

3.1 Geometrical configuration of cooling device

If a larger stiffness coefficient of bearing pedestal is required, particle X ₁ in Table 3 can be used as the final optimization parameter. According to equation (16), the acceptable range of stiffness k_cb of the cooling device is optimized space [1.37 × 10⁵,  2.53 × 10⁵]. The orifice restrictor has been often employed in actual gas bearings, because of its simplicity. The cooling device detailed in Figure 4 is designed for operation at high rotor speed (maximum of 60 krpm) and with impeller temperatures to a maximum of 873 K. The entire cooling device is divided into two hydrostatic zones by the pressure relief groove, which act on the two nodes of the shaft respectively. Each zone has eight equally distributed orifices. In addition, the pressure relief groove with eight inlets connected to the discharge channel to avoid the cross stiffness caused by the dynamic pressure effect and to increase heat dissipation capacity.

In order to prevent the shaft being stuck due to thermal expansion in the shutdown stateand over positioning, the radius clearancebetween the cooling device and the shaft is set to a larger value c_r = 0.05 mm. For lubrication flow, Reynolds number is [16]$$Re=\frac{\rho u{L}_2}{\mu }{\left(\frac{c_r}{L_2}\right)}^2$$(31)where u = 78.54 m/s is the rotor working speed and L₂ = 10 mm is the width of the hydrostatic zones. Consider the extreme case: the density and dynamic viscosity of CO2 near the critical state are ρ_c = 467.6 kg/m³ and μ_c = 2.88 × 10⁻⁵ p_a ⋅ s, respectively. The value of equation (29) is derived as Re = 318.8, which is less than the critical Reynolds number [16] Rc = 1000 of the eccentric ring. Even under extreme operating conditions.

The distance from the geometric center of the rotor to the geometric center of the cooling device is recorded as e_c. In a later section of this paper, the cooling device was designed, and its dynamic characteristics and thermal management were simulated. The carbon dioxide in the radial clearance of hydrostatic zones is still in laminar flow. Table 4 shows the other principal dimensions of SCO2 cooling devices treated in this paper.

Fig. 4 SCO2 cooling devices treated in this paper.

3.2 Orifice diameter of cooling devices

To numerically obtain the static and dynamic characteristics, CFD are employed to predict experimental stiffness coefficients. In this paper, the commercial software Fluent were adopted to calculate the stiffness coefficients. The pressure relief groove in Figure 4 communicates with the atmosphere, therefore the pressure at the hydrostatic zone ends in the radial direction is set to the ambient pressure. The initial pressure boundary conditions at the entrance of the inlet hole is given, and assuming the pressure p_s is 7 MPa.

Figure 5 depicts the recorded dynamic stiffness coefficients versus orifices with diameters of 0.1, 0.15 and 0.2 mm. Different from the conventional fluid lubricated bearing [14], the dynamic pressure effect is minimized by the existence of the pressure relief groove. No significant effect on the stiffness of the cooling device are ascertained between the rotor speeds operating in the range of 0 to 60 krpm. In contrast, the stiffness coefficient of the cooling device changes with the increase of rotor displacement, due to the flow of the orifices is restricted by rotor surface. In the rotor system studied in this paper, the cooling device assumes auxiliary support and heat dissipation, but does not bear the load of the rotor. Note that the bearing capacity does not need to be considered. Therefore, according to optimization space k_cb mentioned above and Figure 6, the diameter of the orifice should be between 0.1 and 0.15 mm. In the following study, the diameters d_s are selected as 0.12 and 0.13 mm.

Fig. 5 Dynamic stiffness coefficient versus rotor displacement with different diameters ds.

Fig. 6 Pressure distributions around an orifice obtained by CFD.

3.3 Dynamic stiffness of cooling device

Figure 6 shows the pressure distributions of the hydrostatic zone. The supply pressure is 7 MPa, and the outlets are set to local atmospheric boundary. Thefigure depicts that the flow of carbon dioxide was mainly restricted at two places where pressure changes rapidly. One is at the entrance of the throttle orifice and the other is at the junction between the throttle orifice and the gap. The gas film at the junction and the throttling effect of the orifice together form the stiffness of the cooling device. Figure 7 shows the variation trend of the static and dynamic stiffness for the cooling device with orifice of 0.12 and 0.13 mm diameter. As indicated in this picture, the stiffness coefficient for rotor speed (n) values of 0 and 50 krpm that are given at the same e_c were almost no change. With the increase of e_c, the stiffness coefficient presents a nonlinear and complex tendency of variation.

Figure 8 shows the relationship between the stiffness coefficient (k_cb) and mass flow (Q) of outlet with orifices of ds = 0.12 mm. The Q here refers to the mass flow rate of the orifice facing the direction of rotor displacement. When e_c is less than 1.0 × 10⁻⁵ mm, the mass flow rate of the orifice remains unchanged, and the stiffness coefficient is also approximately a constant value. Under this condition, the stiffness of the bearing is the inherent rigidity of the CO2 gas film at the junction of the radius clearance and the throttle orifice, which is approximately linear. For 1.0 × 10⁵ m < e_c < 2.5 × 10⁵ m, the stiffness coefficient gradually increases with thegradual decrease of mass flow. For e_c > 2.5 × 10⁵ m, the bearing stiffness curve becomes tortuous because carbon dioxide flow was choked at the outlet of the orifice. As demonstrated in the figure, the mass flow rate of the orifice decreases rapidlyfor this stage.

In order to analyze and verify the influence of the stiffness of the gas bearing on the dynamic characteristics of the optimized rotor system, the inflection points K _c (K _c = [1.46E5,   1.88E5,   3.26E5,  3.46E5,   1.84E5, 8.02E5]) of the stiffness curve in Figure 8 were extracted. Replace the k_ge in X ₁ with K _c to obtain a new particle swarm X _ci (i = length (K _c)). The commercial software ANSYS was employed to establish the finite element model with particle swarm X _ci and experience particle X _e as the parameters. As shown in Figure 9 in the finite element three-dimensional model, the impeller is equivalent to a mass point with the same mass and moment of inertia. The effective vibration mass of the flexible bearing pedestal is m_hub, which is shown as a circular hub in the figure. The inner and outer sides of the hub are virtual oil-lubricated tilted pad bearings and stiffness damping system respectively. The harmonic response analysis results of the finite element model for parameters X _ci and X _e are displayed in Figure 10, corresponding to the solid and dashed lines in the figure, respectively. Compared with the dashed line, the resonance peak corresponding to X _ci has a significant attenuation, and the average attenuation rate is 60.7%. Furthermore, the results of the finite element analysis revealed that the value of K _c was selected in a wider range based on the characteristic number k_ge, rather than limited to the optimization space (Eq. (16)), the resonance peak of the rotor system is still effectively controlled. Which further verifies the effectiveness and robustness of the optimized solution obtained by the CIMPSO algorithm.

Fig. 7 Rotor displacement ec versus stiffness kcb.

Fig. 8 Stiffness coefficient kcb versus mass flow rate Q of the orifice.

Fig. 9 Three-dimensional finite element model of the rotor system.

Fig. 10 Verify the effect of the stiffness coefficient of the gas bearing with orifice for ds = 0.12 mm on the synchronous responses of the rotor.

3.4 Thermal management

Vertical impingements of single or an array of gas jet on the surface can increase the coefficient of convective heat transfer. To this end, for a fixed diameter of 30 mm rotating shaft, a set of different jet-to-target spacings (H_r/d_in = 0.0125 ∼ 5.0) were investigated, and for a fixed diameter of jet inlet (d_in) of 3 mm, three different impeller temperature of 873, 973, and 1073 K were investigated. Figure 2 depicts the current configuration of the cooling device and its structural diagram. There are eight inlet channels with diameter d_in evenly distributed in the circumferential direction of the cooling zone for supplying SCO2. H_r in the figure represents the radius gap of the pressure relief groove. High-density carbon dioxide (CO2) is injected into the pressure relief groove, and then recovered through the discharge channel. Note that the impeller acts as a steady source of thermal energy to heat the rotor-bearing system, and assuming that the temperature of the impeller is heated to 873 K and the supply pressure of CO2 at the inlet channel is 7 MPa. While considering the endothermic effect of gas in the expansion process, Figure 11 depicts the recorded temperature of CO2 versus free end temperature T_c of the rotor. The results reveal a evident reduction in the temperature of the free end relative to that of CO2 at the inlet with d_in = 3 mm and d_in = 4 mm, San Andres et al. obtained similar conclusions with experiments [17]. On the other hand, Figure 11a shows that an inlet channel with a diameter of d_in = 2 mm cannot supply sufficient flow to ensure that the free end temperature T_c equal to or lower than T_CO2.

In SCO2 turbines, an effective thermal management strategy must be ascertained to maintain an acceptable temperature, not affecting significantly the material and lubricating medium properties of the bearings and avoiding axial excessive thermal gradients of the rotor. Therefore, CO2 of higher temperature should be selected as the cooling medium. The cooling device is divided into three zones: the cooling zone in the middle and the hydrostatic zones on both sides. The convective heat transfer occurs mainly in the wall jet area located in the cooling zone. As shown in Figure 12, the overall temperature distribution diagram with T_CO2 = 350 K, the lowest temperature in the cooling zone islocated at the outlet (wall jet area). This indicates that the expansion and endothermic effect of high-density carbon dioxide needs to be carefully considered. Comparing with Figures 12a and b, the temperature gradient, over the whole shaft segment, increases because of the hydrostatic zone. The Figure 12c shows curves denoting the shaft segment temperature gradient with different orifice diameters. The solid line representing the rotor axis temperature in Figure 12b is regarded as a base line since the cooling device does not include any hydrostatic zone. The axial thermal gradient from the heated end toward the free end of the rotor increases due to the hydrostatic zone. The temperatures on free end, on the other hand, does not change significantly with the diameter d_s of the orifice in the hydrostatic zone (Fig. 12).

In terms of geometry, the convective heat transfer coefficient is mainly affected by the jet-to-target spacing (h_r = H_r/d_in), when the gas jet flows into the static environment from a circular section nozzle with a diameter of d_in. In order to investigate the influence of the dimensionless jet radius (H_r/d_in) of the cooling device on the cooling efficiency under the two boundary conditions of a given pressure and a given mass flow rate respectively, the cooling efficiency is evaluated by the temperature of the rotor after cooling. The endothermic effect of gas expansion and the influence of H_r/d_in can be considered in CFD calculation of thermal coupling analysis. It can be concluded from Figure 11 that the inlet with d_in ≥ 3mm can supply enough CO2 to take away the heat transferred from the impeller to the shaft. Therefore, the influence of gas flow on heat dissipation efficiency can be eliminated by setting d_in equal to or greater than 3mm. The heat dissipation efficiency was able to be quantified by the free endtemperature T_c.

Figure 13 concludes with the relative heat transfer performance of different configurations (H_r/d_in = 0.0125 ∼ 5.0) by comparing the average temperature (T_c) of the free end of the rotor. The above temperature parameter was chosen to signify the cooling effect of the cooling configuration for a given heat input (T_s) and nozzle boundary conditions. The inlet temperature of carbon dioxide at the nozzle is 320 K. At the different fixed T_s, Figure 13a depicts the temperature of the rotor after being impacted by the CO2 flow versus radius gap H_r. It is noted that the strong effects of the jet-to-target spacing on the rotor thermal performance at a given mass flow rate of the nozzle. In the interval 0.0125 < H_r/d_in < 0.07 configuration, the same coolant takes away more heat and makes the temperature of the rotor lower than otherwise. As a result, the temperature at free end for that the jet-to-target spacing was varied from 0.1 to 0 case decreases at a faster rate as compared to the larger jet gap (H_r), hence higher efficiency in heat transfer is expected for the much small jet gap configuration. However, the gain in heat transfer comes at a cost of increased pumping power in order to maintain the constant mass flow condition as the other configurations. With the variation of the H_r/d_in from 0.0125 to 5.0, after a fluctuation, T_c continues to increase approximately linearly until reaching the inflection point, and the inflection point is near H_r/d_in = 1. It can be inferred that the heat dissipation efficiency of the cooling device will no longer be significantly attenuated when H_r/d_in > 1, because after the inflection point, the value of T_c will no longer fluctuate significantly with the variation of H_r/d_in.

Globally averaged Temperature at the free end of the rotor has been plotted with respect to associated source temperatures (T_s = 873 K, T_s = 973 K, T_s = 1073K) for all the configurations in the range of 0.0125 ≤ H_r/d_in ≤ 5.0, as shown in Figure 13b. The cooling capacity have been studied for the given pressure of nozzle. Configurations with narrow jet spacings had higher flow resistance due to the presence of viscous force. In the range of 0.0125 < H_r/d_in ≤ 0.17, with increasing H_r/d_in, this CO2 flow resistance reduced, and the temperature at the free end of the rotor approached a linear decreasing trend. For the 0.17 < H_r/d_in ≤ 5, the temperature at the free end fluctuates only slightly, and overall it is close to the horizontal trend, therefore, the best combination of H_r/d_in should be within this range. This paragraph presents a comparative analysis between studies related to jet impingement onto rotor with a given mass flow rate and the studies that deals with jet impingement onto rotor under a given pressure. To this end, the mass flow rate of various jet-to-target spacing H_r/d_in configurations in the study were compared with the temperature T_c of the free end as shown in Figure 13b. Above configurations involving increased mass flow in significant enhancement in heat transfer when compared to the small jet-to-target spacing configuration.

Fig. 11 The temperature of SCO2 versus free end temperature, calculation condition: (a) din = 2 mm, (b) din = 3 mm, (c) din = 4 mm.

Fig. 12 Temperature distribution cloud map and gradient comparison chart, din = 4 mm, Hg = 1.5 mm, TCO2 = 350 K, PCO2 = 7 MPa.(a) With hydrostatic zone. (b) No hydrostatic zone. (c) Temperature gradient comparison.

Fig. 13 Rotor temperature demonstrates that characteristic parameter Hr/din versus cooling efficiency, (a) Given mass flow (q = 0.1 kg/s) of nozzle, (b) Given pressure (1 MPa) of nozzle.

4 Conclusions

Increasingly often in today's world of turbine design and development programs, stakeholders in rotor design desire the use of intelligent optimization algorithms in the design process. A primary reason is that intelligent optimization algorithm have been shown to optimize dynamic performance of rotor, such as reduced vibration and adjustment of critical frequency, more accurately and efficiently than Design experience. In this paper, a multi-objective optimization design is carried out for the rotor system of the SCO2 microturbine. Within this contribution an optimization design procedure of rotor system with emphasis ontwo nonlinear elements is proposed, and the chaotic interval particle swarm optimization (CIMPSO) algorithm is summarized.

This algorithm is used to optimize the rotor system of SCO2 turbine for multiple times, and the global non-inferior optimal or suboptimal solution sets are obtained. Commercial software ANSYS was used to establish a finite element model base on the solution sets for simulation verification. The results show that the resonance peaks of the rotor are effectively attenuated, and the critical speed is far away from the operating speed point, which verifies the robustness of the algorithm.

Demonstrated cooling device at high temperature is fundamental to enable the implementation of active magnetic bearing (AMB) and Oil-lubricated tilted pad bearings (OTB) into SCO2 turbines applications. The cooling device is simulated and the results revealed thatSCO2 cooling flow rates of high pressure character ensure an adequate thermal management for turbines operation at high temperatures. Determination of the high cooling efficiency for adequate thermal management is an issue of current research. For the coolant temperatures, six different jet-to-target spacing (H_r/d_in = 1/8, 1/6, 1/4, 1/3, 1/2, 3/4) were investigated, and for the fixed coolant temperature (T_CO2 = 320K), two different nozzle boundary conditions of given inlet mass flow rate (q = 0.1 kg/s) and given inlet pressure (p_s = 1 MPα) were investigated. Flow and thermal field computations are performed for the above two boundary conditions to investigate the effects of jet-to-target spacing (0.0125 ≤ H_r/d_in ≤ 5.0) and other parameters on the heat transfer from the impingement surface. The following conclusions about the heat transfer were drawn from the study.

• For high-pressure and high-density carbon dioxide jet impinges heat transfer, the jet mass flux and the heat absorption by gas expansion have significant effects on the heat transfer characteristics. Jet impingement heat transfer for given the boundary conditions of the pressure and temperature of the CO2 coolant was primarily affected by the mass flow rate, while the jet-to-target spacing H_r/d_in effect was secondary. In the range of 0.0125 ≤ H_r/d_in ≤ 0.33, with larger H_r/d_in yielding in higher heat transfer. This trend was observed in a variety of source temperature (T_s = 873 K, T_s = 973 K, T_s = 1073 K) boundary conditions.

• When heat transfer at constant pumping power (given pressure boundary conditions) was considered, the 0.33 ≤ H_r/d_in case was the best performing configuration, while H_r/d_in < 1 performed the best under the given mass flow rate conditions.

• For jet impingement heat transfer with a given mass flow rate, the relationship between heat transfer efficiency and jet-to-target spacing is a highly nonlinear function. At low jet-to-target separation distance (0.0125 < H_r/d_in < 0.1), the global heat transfer characteristic decrease almost linearly with increasing jet-to-target spacing, while at intermediate separation distance (0.1 < H_r/d_in < 0.17), the global heat transfer capacity initially increases mildly and then attenuates approximately linearly (0.17 < H_r/d_in < 1.0) with increasing H_r/d_in.

Funding

The research has been funded by National Key Projects of China (No. 2018YFB2000501).

The authors gratefully acknowledge the financial support providedby the National Natural Science Foundation of China (No. 51375148).

The research has also been funded by Henan Provincial Wisdom Introduction Project (Outstanding Foreign Scientist Studio) (No. GZS2018005).

Acknowledgments

The authors would like to acknowledge the support fromHenan 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

The ZDT problem proposed by Zitzler is a commonly used two-objective test function [18]. In this paper, ZDT1 and ZDT2 are selected as the tested objects:$$\begin{array}{l}{f}_1\left(x\right)=x_1\\ f_2\left(x\right)=g\left(x\right)\left[1-{\left(f_1\left(x\right)/g\left(x\right)\right)}^m\right]\\ g\left(x\right)=1+9\left({\displaystyle \sum _{i=2}^n}{x}_i\right)/\left(n-1\right)\\ x_i\in \left[\mathrm{0,1}\right],i=\mathrm{1,2},\cdots ,n\\ n=30\end{array}$$(A.1)

The parameters m = 1/2 and m = 2 correspond to ZDT1 (convex function) and ZDT2 (concave function) respectively. Their ideal Pareto front has been given, as shown by the dotted and dashed lines in Figure A.1. The parameter x₁ affects both sub-objective functions and is therefore selected as a characteristic number element. The blue line segments in the figure are the solution intervals obtained by numerical experiments. The number of initial particles corresponding to each numerical experiment is 50, the number of iteration steps is 200, and the total number of experiments is 200. Compared with the ideal Pareto front, it is obvious that the experimental solution converges to the non-inferior optimal solution (PF true) with a certain error range. The average errors of the test solutions of ZDT1 and ZDT2 are e₁ = 16.1% and e₂ = 11.69%, respectively.

Fig. A.1 Ideal Pareto front (PFtrue) and test solution set (EXPERIMENT).

References

All Tables

All Figures

	Fig. 1 Supercritical CO2 radial-in-flow turbine lay-out.
In the text

	Fig. 2 Approximate rotor shaft geometry.
In the text

	Fig. 3 Rotor imbalance response before and after optimization design.
In the text

	Fig. 4 SCO2 cooling devices treated in this paper.
In the text

	Fig. 5 Dynamic stiffness coefficient versus rotor displacement with different diameters ds.
In the text

	Fig. 6 Pressure distributions around an orifice obtained by CFD.
In the text

	Fig. 7 Rotor displacement ec versus stiffness kcb.
In the text

	Fig. 8 Stiffness coefficient kcb versus mass flow rate Q of the orifice.
In the text

	Fig. 9 Three-dimensional finite element model of the rotor system.
In the text

	Fig. 10 Verify the effect of the stiffness coefficient of the gas bearing with orifice for ds = 0.12 mm on the synchronous responses of the rotor.
In the text

	Fig. 11 The temperature of SCO2 versus free end temperature, calculation condition: (a) din = 2 mm, (b) din = 3 mm, (c) din = 4 mm.
In the text

	Fig. 12 Temperature distribution cloud map and gradient comparison chart, din = 4 mm, Hg = 1.5 mm, TCO2 = 350 K, PCO2 = 7 MPa.(a) With hydrostatic zone. (b) No hydrostatic zone. (c) Temperature gradient comparison.
In the text

	Fig. 13 Rotor temperature demonstrates that characteristic parameter Hr/din versus cooling efficiency, (a) Given mass flow (q = 0.1 kg/s) of nozzle, (b) Given pressure (1 MPa) of nozzle.
In the text

	Fig. A.1 Ideal Pareto front (PFtrue) and test solution set (EXPERIMENT).
In the text