Using an optimisation strategy to design a supercritical CO2 radial inflow turbine transonic stator

During the operation of supercritical CO (sCO ) radial inflow turbines (RITs), sonic conditions may happen, which will decrease their efficiency. Non-standard design geometries are needed to reduce the losses. However, modifying the turbine geometry is a challenge, especially as the optimum shape may be non-intuitive; commercial computational fluid dynamics (CFD) software is not friendly towards automation by surrogate scripts, and the traditional design method for blade passage shape is not suited to the optimisation of aerodynamic turbine design. Hence, in this study, a stator nozzle for a 120 kW sCO RIT, whose geometries are obtained from in-house preliminary design code TOPGEN , has been optimised with a modularised geometry optimiser. To parametrise the geometry, a parametrised stator mesh generator based on a meridional streamline is developed. Bulk CFD simulations are carried out with OpenFOAM , to form a Pareto front with a vector of 14 variables. Three final stators have been selected: they are the optimised Mach number distribution ( ), the optimised outlet flow angle distribution ( ) and a compromise case. The the outlet boundary properties of these stators are extracted and discussed. The optimised case has the shortest divergent nozzle, returns the best and the least loss; the optimised case has the longest divergent nozzle and returns a better-uniformed outlet flow; the compromise case has the medium length of the divergent case. All these stators have fixed mass flow rates to allow de-coupling of the upstream system from the rotor. This study will ultimately benefit the development of the sCO power cycle.


Introduction
With the rapid development of human society, a large quantity of pollutants and global warming gas have been released into the atmosphere, which is bringing about environmental crisis. This has lead to a firm determination to save energy and reduce emissions. The application of renewable energy will be the ultimate solution for people to reach these goals. There are several renewable energy sources available, such as biomass, wind energy, hydraulic energy, nuclear energy and solar energy. Among these renewable energy sources, solar energy is the most abundant and it has already been used by human society in various forms since ancient times. However, today solar energy provides less than 0.53 % of the total global energy (Mekhilef et al., 2011). This creates an opportunity to develop solar energy technologies.
CONTACT Xujiang Wang x.wang@sdu.edu.cn Currently, a significant focus has been placed on the development of solar thermal energy applications. One of the most efficient and economic paths is to use concentrated solar power (CSP) systems. Here, solar energy is collected and converted into heat using a field of heliostats and a thermal receiver, and furthermore the use of cost effective thermal energy storage can be considered (Harries et al., 2012). Next, a thermal heat engine (power cycle) is used to convert thermal energy into electric energy.
Owing to several advantages of CO 2 -for example it can easily reach the critical point (7.38 MPa,304.25 K) compared to H 2 O (22.064 MPa, 647.1 K) or other fluids (Dostal, 2004) and rapid changes in density reduce compressor work and lead to higher heat transfer coefficients in the recuperator and pre-cooler -the supercritical carbon dioxide (sCO 2 ) power cycle is considered as a promising power cycle for the next generation. As proposed by Dostal (2004) the Brayton sCO 2 cycle can offer higher efficiency than the Rankine steam cycle when the turbine inlet temperature is higher than 550 • C.
Turbomachinery is the most important part of the sCO 2 power cycle, whose performance affects the overall cycle efficiency significantly. Brun et al. (2017) have shown that about every 2% increase in turbine efficiency results in approximately 1% increase in cycle efficiency, while the influence of compressor efficiency is approximately half that. sCO 2 axial turbine development is aimed at future large-scale energy production (>50 MW), where sCO 2 can replace steam or other working fluids in the mid to long term (Southwest Research Institute, 2019). However, sCO 2 also has potential at smaller scales. For power cycles in the range 0.1 to 25 MW, using sCO 2 allows the use of more efficient radial inflow turbines (RITs). Using an RIT architecture has several advantages, such as fewer seals, rotors that are more vibration resistant and better suited for the high blade loadings, and simpler overall system architectures to name a few (Qi et al., 2017;Ventura et al., 2012). This is the case owing to the combined effects of the highly dense working fluid and comparably low flow rates, resulting in highly power-dense machines. In this study we will focus on RIT designs.
To reach higher efficiency, a good design methodology is required for sCO 2 RITs. Currently, sCO 2 RITs or turbine components are designed following traditional procedures, i.e. using quasi-one-dimensional inviscid solvers combined with time-consuming computational fluid dynamics (CFD) simulations, which apply the thermodynamics models available at the time (Pasquale et al., 2013;Qi et al., 2017). Relying only on these designs will lead to risks for the RITs, especially for RITs operating on a high-pressure ratio, for which supersonic loss will be detrimental to the overall efficiency. To maximise the efficiency, one possible option is to use non-standard geometries that optimise blade geometry to suit the flow inside the blade channel. Owing to the complexities of the flow and the 3D nature of the passages, modifying the shape is challenging, especially as the optimum shape may be non-intuitive.
Multiple tries to find the optimum reduce the possibility of experimental research. An alternative approach intensively uses CFD with appropriate optimisation strategies, which allows for more extensive exploration of the geometry design spaces. As reported in our previous study (Qi, Xu, Han, & Zhang, 2022) ('Development and application of a modularized geometry optimizer for future supercritical CO 2 turbomachinery optimization'), a modularised geometry optimiser based on the Nelder-Mead algorithm was developed and validated.
As bulk CFD simulation needs to be automatically updated and evaluated during the optimisation progress to explore the geometry design spaces, commercial CFD software is not usually automation friendly regarding surrogate scripts. Hence, using open-source CFD software is a good option. OpenFOAM is a leading open-source project for continuum mechanics and CFD applications that has some support for non-ideal computational fluid dynamics simulations and a flexible framework for future developments (Weller et al., 1998). In this study, this optimiser will be applied with the open-source CFD software OpenFOAM to carry out the optimisation.
To optimise the blade shape, a strategy for geometry parametrisation is needed (Da et al., 2020;Sui et al., 2021). Traditionally, the shape of blade passages on an RIT is defined using a hub contour, a shroud contour, a wrap angle and a local blade angle β, whose principles of definition are all about the passage meridional length. Such a geometry definition is logical methodology from a structural, manufacturing and mechanical design point of view, but is less well suited for the optimisation of aerodynamic turbine design (Jahn & Jacobs, 2016). The study of Bonaiuti and Zangeneh (2009) has shown that using inverse design methods for geometry definition leads to reduced complexity in correlations between design parameters and performance parameters, which then allows more efficient analysis and optimisation. For the aerodynamic optimisation of RITs, the ideal design parameters are local flow direction and rate of change of flow area (the area normal to the flow direction). But in fact, the underlying physical changes, such as the differences in the rate of turning angle and passage length, affect the blade loading and profile losses, resulting in a more gradual expansion and possibly more subtle changes to the flow field in the blade channel. Thus, to gain an understanding of the dominant flow phenomena affecting performance for optimisation, it is most effective to parametrise the design space using parameters that directly relate to flow performance and operation. In this study, a parametrised mesh generator that can easily be connected and controlled by the optimiser is also developed.
Considering a whole turbine stage, the stator plays an important role. For a small sCO 2 RIT operating at a high rotational speed, it is desirable to keep the mass flow rate (ṁ) constant and steady in order to maintain the requirements of the closed loop cycle. If applying a convergent nozzle, fluctuation of the downstream stator exit or rotor inlet conditions will highly affect the wholeṁ cycle. An unstableṁ will highly affect the stability of the whole cycle. There are several methods of maintaining a constantṁ for RITs. One way to achieve a constantṁ is to use a transonic stator nozzle, as once the flow reaches a chocked condition at the nozzle throat, theṁ becomes fixed.
However, expansions through the sCO 2 RITs stators are commonly characterised by large pressure ratios, and a low speed of sound (a) of the working fluid is generally observed. The flow can easily reach sonic conditions, and this increases viscous (friction) losses in the stator. In this case, non-standard geometries such as stator blade shapes with converging-diverging channels are therefore required to control theṁ correctly and to provide a better rotor inlet flow field. In contrast with optimising the whole turbine stage, which is far more complex today, this study is focused on the optimisation of the stator nozzle of a sCO 2 RIT.
Before showing our research objectives, here list of the major difficulties and challenges encountered during the research: (1) modifying the turbine blade shape is challenging especially as the optimum shape may be nonintuitive; (2) commercial CFD software is not friendly with automation by surrogate scripts; (3) traditionally defined method of blade passage shape is not suited for the optimisation of the aerodynamic turbine design.
Based on the above discussion, this study will carry out optimisation for the sCO 2 RIT stator, and the objectives for this study are: (1) to develop a parametrised turbine stator mesh generator, which could generate the stator mesh based on different stator geometric parameters; (2) to optimise the turbine stator base-line geometry towards a different combination of weighting factors on objective function; (3) to use open-source CFD software OpenFOAM to acquire the flow field within the blade channel of the stators; (4) to discuss the aerodynamic performance of the optimised turbine stator geometries.

Motivation for optimisation
In previous study (Qi et al., 2017), preliminary designs for 100 kW to 200 kW sCO 2 RITs are developed with in-house code TOPGEN (Ventura et al., 2012). It would be well suited if the target power is within this range. Thereafter a slightly large output power 120 kW turbine is preliminarily designed by TOPGEN, and the outlet flow variables for the stator are given as the targets for the optimisation problem. Following that, a non-standard designed transonic nozzle, which can return both correct mass flow rate (ṁ) and require velocity triangle will be designed and discussed. The detailed geometry parameters and targets for initialising the optimiser are provided in Table 1. Hence, in this paper, the optimisation process for the 120 kW turbine transonic stator will be presented. In the meanwhile, the optimised stator geometries will be discussed and analyzed.

Methodology
In this section, the optimisation methodology for a convergent-divergent stator nozzle is provided.

Objective function
It is essential correctly to define the objective function for an optimisation problem, which delivers the targets to the optimiser. Considering that the outlet flow fields are the most important for the stator, as they affect the operation of the downstream rotor, the objective functions are set based on the stator outlet flow properties. Figure 1 presents an example of the stator geometry and flow fields. In this project, six different targets (highlighted in bold font in Table 1) leading to a multi-objective optimisation problem need to be reached, and the linear combination of objects can be addressed as where should be minimised. ϕ Vi is the cost term related to different targets, and W i is the corresponding weight.
Obviously, the total mass flow rateṁ 0 is one of the cost terms, which is related to the turbine output power and cycle operating conditions. The mass flow rate for a single turbine blade passage isṁ =ṁ 0 /N. Stator efficiency needs to be considered, and is defined as the ratio of total pressures p 0 , i.e. p 01 /p 00 . To achieve the correct spouting velocity, the magnitude of the velocity and the flow averaged angle are critical. Considering the local sound speed a, the average Mach number M is selected as one target. The average outlet flow angle α is also considered. Besides M and α, their deviations are also important, as they reflect the stability of operation. Therefore, the standard deviation of the Mach number, σ M , and outlet flow angle, σ α , are added as targets. In all, there are six objective terms that need to be considered, resulting in the total objective function as The standard deviation of a variable (x) is calculated through the following equation: where N is the total number of samples and x is the mean value of variable x. The standard deviation gives a measurement of the deviation from the expectation or average value. Among the six objective terms, M, σ M , α and σ α can be modified into two equations, as ϕ M σ and ϕ α σ . They correspond to the deviations from the target values rather than the expectation or mean value: Here, be careful to note that M σ is different from σ M . The former indicates the deviation of the Mach number from the target, and the latter indicates the standard deviation of the Mach number distribution on the outlet patch; likewise α σ and σ α . This reduces the number of objective terms to four, which simplifies the complexity significantly. The detailed equations are listed as follows.
• ϕ M σ is the deviation between the flux averaged Mach number and the target Mach number (marked with the subscript 'tar'), calculated by • ϕ α σ is the deviation between the flux averaged outlet flow angle and the target outlet flow angle, • ϕṁ is the term related to the mass flow rate, which is calculated by • ϕ p 0 is the total pressure cost term, which is calculated by where ρ i is the fluid density, u mi is the meridional component of the flow velocity, A i is the area, M i is the Mach number of the outlet boundary patch and n i is the normal vector of the face. n is the total face number of the outlet boundary patch. The cost contributions of Equation (4) are evaluated at the outlet boundary of the computational domain. The modularised optimiser will minimise the value of the objective function. One thing to be considered is that, if looking at the objective terms, it can be seen that all the terms are odd functions apart fromṁ. The objective function forṁ is a parabolic function, and that means it has a lowest point. Hence, for a given optimisation problem, a large weighting factor forṁ is given first to force the optimiser to find the correctṁ. Then the cases with the correctṁ will be used as an initial step to start the simplex optimisation process for other objective terms.

Stator geometry generation strategy
In the process of optimising the turbine stator blade geometry, the most important part is geometry parametrisation, as poor parametrisation can result in too many degrees of freedom or geometries that are not flexible. In this study, a 2D schematic of the stator geometry with multiple control points, curves and blocks is shown in Figure 2. The stator geometry is generated through the following procedures.

First: Set the fixed parameters
The basic parameters for the transonic stator nozzle are provided in Table 2, which fix the overall dimensions for the studied nozzle. These parameters are obtained from TOPGEN (Ventura et al., 2012). It should be noted that TOPGEN does not account for the stator geometry design (except that the stator outlet radius R 2 and the outlet flow angle α 2 are fixed); thus, this is only used as a starting point.
By setting the inlet and exit radii, the nozzle ring domain is set. The outlet radius determines the blade outlet radius, and the center of the trailing edge. The vane number sets the pitch angle of the single passage, as θ p = 2π/Z r . The stator start angle θ s determines the starting angle relative to zero radians in an anti-clockwise direction of the blade pitch. With these five parameters, the position of the stator blade can be set.

Secondly: Draw the floating control points
Once the three parameters are determined, the other points to control the blade shape need to be determined. As shown in Figure 2, the floating control points are drawn as red dots. For example, the points XR2 and XR1 control the positions of the two circles (as shown in green, which control the nozzle throat wall curvature), and the point Xt controls the position of the nozzle throat section.

Thirdly: Generate the curves
The nozzle shape is composed by Bézier curves, as shown in Figure A1 of Appendix 1. Three different Bézier curves are generated using floating geometry control points. The control points for drawing the Bézier curves are marked with different colors and '1', '2' or '3' as subscripts. Line 1 and line 2 control the outer part of the nozzle, while line 3 controls the other parts of the blade.

Mesh blocking scheme
The stator represents a most challenging station in terms of high-quality mesh generation due to its geometric characteristics, namely high outlet angle, small passage areas (between vanes in the pitch-wise direction) and small radius at the trailing edge. This section presents the mesh blocking structure for transonic RIT stators. The geometry and control lines are drawn for the stator as discussed in Section 3.2. To save computational cost and as the stator flow is primarily 2D, a 2D mesh is used. The stator mesh is committed to a 2D simulation, and thus the mesh details in the third dimension can be omitted in this project. As shown in Figure 2, the floating points to control the geometry are marked with black dots. The fluid domain is decomposed into 30 blocks, and the block names are marked in blue. After adjusting the cell density and the grid gradient, an example of the resulting mesh derived from this blocking structure is presented in Figure 3. A detailed stator geometry parametrisation strategy will be discussed in Appendix 1.
In Figure 3, the boundary patch is defined as 'Inlet', 'Outlet', 'Inside wall' and 'Outside wall'. Flow comes into the stator nozzle from the 'inlet' patch and accelerates out from the 'Outlet' patch. Rectangular cells are distributed and filled from the throat to the divergent part of the nozzle. Near-wall clustering is applied to the mesh with the Roberts cluster function (Jacobs et al., 2015;Roberts, 1997) to keep an optimum y + value, which allows proper wall functions to be used. A high quality mesh around the trailing edge ensures correct simulation of the wakes, as the separation from the trailing edge may affect the downstream flow field significantly.
As the optimisation of the geometry needs iterations from mesh generation, CFD simulation, results evaluation and again next iteration from mesh generation, thus the automation of mesh generation is important. The mesh is generated with the mesh tool package of Eilmer3 (Jacobs et al., 2015) since the mesh generator of Eilmer3 is well-developed, open-source and allows easy access for Python script. For an optimisation problem, easy access to the mesh generator by control script is essential, which allows quick updating of the mesh. With the help of the e3prepToFoam (Jahn & Qin, 2015) utility, the generated mesh is converted to OpenFOAM format. The mesh reports indicate that this stator mesh has a fairly good quality: the maximum aspect ratio (875.97), maximum non-orthogonality (69.5) and maximum skewness (3.19) are all below the threshold value of OpenFOAM.

CFD simulation details
The stator boundaries are depicted in Figure 3. These boundaries are labeled as: Inlet (OF_inlet_00), Outlet (OF_outlet_00), outside wall (OF_wall_00), inside wall (OF_wall_01) for OpenFOAM. The corresponding boundary conditions are presented in Table 3.
In this project, the inlet boundary is taken from the test case and the outlet boundary is determined by the leading edge of the following rotor blade. For the present situation, the periodic boundaries of the passage are  separated into two parts: one is located from the inlet to the leading edge and the other is located from the trailing edge to the stator outlet. The periodic boundaries are marked as OF_cyclic_00 and shown as blue lines in Figure 2. The inlet total temperature is 833.15 K and the inlet total pressure is 20.0 MPa, which is calculated from the power cycle requirements (Qi et al., 2017). The stator outlet pressure is 13.0 MPa, which leads to a pressure ratio of 1.538. The solver used to carry out the simulation is RGDFoam, which was developed and validated by Qi, Xu, Han, and Jahn (2022). This solver was enhanced with the ability to solve non-ideal computational fluid dynamics (NICFD) fields from an approximate Riemann solver (HLLC) by Borm et al. (2011). To validate the solver in NICFD fields, three cases are simulated (Qi, Xu, Han, & Jahn, 2022). First, a 2D simulation of a NASA transonic air nozzle operating with air is used to confirm the solver's ability to simulate correctly transonic flow phenomena and shock waves. Then, the simulation of a 2D VKI cascade operated with MDM gas in a stationary frame is used to assess the solver's ability to simulate correctly non-ideal gas flows typical in industrial applications. Lastly, the simulation of dense gas flow (MD 4 M) passing a backward ramp to illustrate the ability of the approximate Riemann flux calculator and look-up table mechanism to work well in the non-ideal region of fluid properties. The concise validation process is presented in Appendix 2, and full a validation process has already been published in a study by Qi, Xu, Han, and Jahn (2022).
The k-ω-SST turbulence model is applied to close the momentum equation. Attention is focused on the first layer mesh from the wall to guarantee that the y + value is less than 30, and to allow wall functions to be employed to capture the near wall and to reduce the computational load. The Runge-Kutta method is applied to accelerate the convergence speed.
To carry out the grid dependency study, six meshes with sizes of 2.1k, 4.1k, 8.5k, 15.6k, 35.9k and 68.2k cells are used. Two approaches are used to compare the flow structures of interest, one being the flux properties along the meridional line of the blade channel and the another being the properties on the outlet patch. Figure 4 shows the schematic for extracting the fluid properties along the meridional line of the channel. The meridional line is defined as the half-way point (in the tangential direction) between the different boundaries, which indicates the flow direction. In this schematic, the meridional line is shown as a black solid line. Ten points are selected along the line, which are used as the base points for the slices. The slice surfaces are then created from the base points and a vector perpendicular to the radial direction. Flux average fluid properties (φ F , which can be replaced by pressure p, density ρ, temperature T, etc.) are calculated for every single surface by where u is the fluid velocity, A is the corresponding surface area, ρ is the density and n is the face normal vector. The subscript F denotes that the fluid properties are averaged using the flux average method. Figure 4 shows the M contours on different slices. This procedure allows evaluation of the φ along the stator passage. Figure 5 shows a comparison of the flux averaged slice M for the six different cases. The M points are extracted from 50 different slices along the meridional line. It can be seen that the difference between the 15.6k and 35.9k case is 2%, while between the 35.9k and the 68.2k case is 0.7%, respectively. The comparison illustrates that the 15.6k case will be a future potential case to carry out bulk CFD simulations under acceptable accuracy.
As the optimiser evaluates the flow properties on the outlet patch to determine whether the design reaches the target or not, the dependency studies should also be carried out based on the boundary patch properties, for example the outlet flow angle, outlet Mach number variation, etc. Figure 6 compares the Mach number variation and outlet flow angle (α) variation across the stator outlet between different computational meshes. It can be seen from Figure 6 that the 35.9k grid returns results whose error is acceptable (the error for the flux averaged Mach number is 4.6%, and for the flux averaged α is less than 1%). A discontinuity is observed in the M and α distributions for all the 8.5k, 4.1k and 2.1k cell grids. This is due to the junction point of three blocks, namely UT, LE and LT2 shown in Figure 2, where discontinuities of face areas appear for sparse cells. To reduce the discontinuity on face areas, increasing the grid density is the only solution. Based on the fact that the 15.6k mesh appropriately resolves the meridional Mach number and the flow features on the outlet patch, the 15.6k cell mesh is used for the following optimisation process.
To make an independent study for the CO 2 density property tables, we first use the tabular data error estimation tool provided in the study by Qi, Xu, Han, and Jahn (2022). This tool was developed to preliminarily estimate the tabular data accuracy in given primary parameter ranges, i.e. p or T. In this study, the range of p and T should cover the potential simulation spaces, hence p uses 10 to 30 MPa, and T uses 700 to 1200 K. The recommended tabular data resolution is 35(T)×15(T), i.e. in the given T and p range, there are 35 T nodes and 15 P nodes to generate the thermophysical and transport properties tables. However, to study the independence of the table properties resolution, four series of look-up tables are generated and used to carry out the simulations. Figure 7 shows the result of the independence study for the property tables. It can be seen that the recommended resolution, 35(T)×15(p), can return acceptable accuracy. Once the resolution is larger than 35(T)×15(p), the simulated flux averaged Mach number will not change. Hence, in this research, a 40(T)×20(p) resolution is used to balance the computational cost and accuracy.

Results and discussion
In this section, the results of several optimised cases are presented.
First, the optimiser is operated with the initial conditions presented in Table 1. A largeṁ weighting factor is set on the optimiser, to find the optimumṁ value. The resulting design is identified as point A in Figure 8. The largeṁ weighting forces the optimiser to find a geometry with the correctṁ, which is labeled as point B in Figure 8, and will be used as the baseline case in the following context. The corresponding geometry is listed in Table 4.
Once the correctedṁ (with ±5% mismatch) is achieved, a new starting simplex will be created to initialise the optimisation of the other parameters. The     following sections present details of the subsequent optimisation.

The optimisation paths
Owing to the nature of the multi-objective optimisation problem, a Pareto front is required to determine the trade-off between different optimisation targets. To reach different local minima, multiple combinations of weighting factors for M σ and α σ are set. The resulting Pareto front is shown in Figure 8. It's worth noting that this Pareto front is unlike the normal one, owing to the limits on the computational resources, and we only have finite tries on different combinations of weighting factors. The dots on this Pareto front figure are CFD results, with color to show the flux averaged total pressure p 0 at the outlet boundary patch. The initial simplex is located in the center of the figure, around point B. Applying different combinations of weighting factors leads the optimiser to walk in different directions. One route is walking in the top left direction, which means the optimiser tends to find a case with lower M σ and higher p 0 , but a loss in performance with respect to α σ . On the contrary, another route for the optimiser is to walk in the bottom right direction, which implies that the optimiser tends to find a case with lower α σ , but with higher loss (lower p 0 ) and higher M σ .
To get a more uniform flow, i.e. to reduce α σ , it requires an increase in nozzle length. However, reducing the M σ requires a decrease in nozzle length, and reducing the friction loss, thus the nozzle reduces the loss in final kinetic energy, which is coincident with higher p 0 values.
To show the optimisation results, three different stator geometries are picked as possibe optimal geometries, marked with a, b and c, as shown in Figure 8. Case a has the lowest M σ value (marked as 'M σ optimised'), Case b has the lowest α σ value (marked as 'α σ optimised') and Case c is a compromise for both M σ and α σ value (marked as 'balanced performance case'). The values of the geometry state vectors for cases a, b and c are shown in Table 4. The geometry profiles, the flow fields and the outlet flow properties are compared in detail and discussed in the following sections. Figure 9 shows a shape comparison between the baseline nozzle and the M σ optimised nozzle (case a from Figure 8). To have a better comparison, the two nozzle sketches are re-centered based on the mid-point of the nozzle throat. The outlet boundary of the calculation domain is kept, to show the original position of the nozzle throat referenced to the outlet. It can be seen from Figure 9 that the optimised M σ nozzle has a shorter divergent section. The nozzle centerline of M σ has a larger angle than the baseline case. To have a detailed view of the flow field, Figure 10 shows the Mach number contours and the normalised total pressure (p 0 /p 01 ) contours for both the baseline case and the optimised M α case.

Optimisation focus on the Mach number deviation, M σ
It can be seen from Figure 10(a) that the flow accelerates from a low Mach number to a supersonic state at the nozzle throat. A maximum Mach number of approximately 1.3 is reached. The streamlines show the flow directions in the nozzle. It can be noticed that a large separation occurs near the outside wall, and the streamlines of the separation show that a large vortex is formed near the outside wall. This is due to the limitedṁ from the given inlet condition, which cannot fully occupy the downstream nozzle. This large vortex pushes the mean stream towards the inside wall, which confirms again that the optimiser can find a correct nozzle shape to constrain theṁ via 'smartly manipulating' the vortex. Figure 10(b) shows the total pressure contours for the baseline nozzle. It can be seen from the figure that the choke flow at the throat generates the largest losses. When comparing the M σ optimised case to the baseline case, it can be seen from Figure 10(c) that the separation region is reduced owing to a shorter divergent nozzle. More uniform outflow is achieved. What's more, the total pressure contours shown in Figure 10(d) show that the total pressure downstream is enhanced, especially for the main flow stream, which is because the shorter divergent part of the nozzle reduces the viscous (friction) losses. Thus, this nozzle has less loss and a higher exit outlet Mach number than the baseline case. Figure 11 shows the distribution of the Mach number, outlet flow angle α and total pressure p 0 /p 01 on the outlet patch, as a function of pitch angle. In Figure 11(a), the M σ optimised case has an average Mach number of 0.62, which is a 4.6% deviation to the target Mach number, 0.65. This is better than the baseline case, whose average Mach number is 0.53, a 18.5% deviation from 0.65. Compared to the baseline case, less deviation occurs for the M σ optimised case. Similarly, Figure 11(b) shows a comparison for the outlet flow angle α distribution. Considering the mean flow angle, the M σ optimised case performs worse than the baseline case; however, the deviation of α is reduced slightly. Finally, Figure 11(c) illustrates the total pressure (p 0 /p 01 ) distribution along the stator outlet. Higher total pressures are observed across all angular positions for the M σ optimised case compared to the baseline case. This confirms that the M σ optimised case has fewer losses than the baseline case.

Optimisation focus on the deviation of outlet flow angle, α σ
Another selected case focuses on the outlet flow angle (α) distribution. It corresponds to case b in Figure 8 and is labeled 'α σ optimised'. Very similar to the nozzle shown in the previous study (Qi, Xu, Han, & Zhang, 2022), the α σ optimised case has a long divergent part, as marked in red in Figure 12. The optimiser stops due to the nozzle length reaching the limit. This means that, to increase the uniformity of the flow, the nozzle length tends to be longer.
The Mach number contours and the normalised total pressure contours for the α σ optimised case are shown in Figure 13. It can be seen from Figure 13(a) that a more uniform flow compared to Figure 13(b) is achieved.
The mainstream attaches to the inside wall of the nozzle. Compared with the Mach contours of the baseline case shown in Figure 10(a), the separation on the inside wall is more or less eliminated. However, the separation region along the outside wall is larger than in the baseline case. Figure 13(b) illustrates the normalised total pressure (p 0 /p 01 ) contour for the α σ optimised case. Compared with the p 0 /p 01 contour for the baseline case, the mainstream p 0 /p 01 at the outlet section of the α σ is higher than in the baseline case. This phenomenon confirms that the optimiser can successfully find the geometry with lower dynamic loss.
To give a better illustration of the outlet patch, Figure 14(a) compares the distribution of M between the α σ optimised case and the baseline case. It can be  seen that the outlet M distribution of the α σ optimised case is slightly better than the baseline case. But the average Mach number (M) is still 13.8% away from the target M tar . Figure 14(b) illustrates that the α σ optimised case has a better distribution of the outlet flow angle α than the baseline case, with only a 2.7% deviation to the target outlet angle α tar , and the deviation is slightly better than in the baseline case. The performance is enhanced as p 0 /p 01 increases, as shown in Figure 14(c).

Optimisation for a balanced performance nozzle
The third case is a balance between optimised M σ and α σ , which is labelled case c in Figure 8. The comparison of the geometry to the baseline case is shown in Figure 15. This nozzle has a slightly shorter nozzle than the baseline case and the outside wall of the nozzle expands a little. Figure 16 shows the contours of M and p 0 /p 01 for the balanced performance case. Compared with the baseline case shown in Figure 10, the separations on both sides of  the divergent part of the nozzle are reduced. The separation near the outside wall is larger than the M σ optimised case, but smaller than the α σ optimised case. The p 0 /p 01 is shown in Figure 16(b), which illustrates that the loss is reduced compared to the baseline case. Compared to the other two optimised cases, the loss for the balanced performance case is between them, which is in agreement with Figure 8. Figure 17 shows a comparison of the Mach number, α and p 0 /p 01 distributions on the outlet patch for the balanced performance nozzle. It can be seen from Figure 17(a) that the average Mach number (M) of the balanced performance case is closer than the baseline case, but the deviation is larger than the M σ optimised case. Figure 17(b) shows the balanced performance case has an enhancement in the mean α and a slight enhancement of the deviation for the distribution, but the mean α deviation is still worse than the α σ optimised case. Similarly, Figure 17(c) shows the p 0 /p 01 distribution for the balanced performance case. Enhancement of the performance can be noticed from this figure. Figure 18 shows a comparison of the three optimised stator nozzle geometries. For easier comparison, all geometries are stacked together according to the position of the midpoint.

Comparison of three different stators
Obviously, the M σ optimised nozzle has the shortest length among these nozzles, and the α σ optimised has the longest divergent part. However, the α σ optimised nozzle has a thinner trailing edge, which may potentially affect the structural strength.
The separation characteristics are different for these three cases. There is a large separation along the outside wall for both of them. The optimised M σ stator has the smallest separation while the α σ optimised stator has the largest separation. It can be seen that small  separations occur near the inside wall for both baseline cases, M σ and balanced performance stators. The α σ optimised case more or less eliminates the separation near the inside wall, and a uniform exit flow is exhibited. The large separations on both optimised cases indicate that the optimiser can adjust the internal fluid field by using a vortex. The geometries are directly manipulated with the optimiser, and the outlet flow fields are perceived by the optimiser. But, limited by the mesh generator and the optimisation algorithm (Nelder-Mead), currently, the optimiser can only find the local minimum for an optimisation problem.

Conclusions
In this study, the optimisation process with multiple techniques for a sCO 2 radial inflow turbine stator is illustrated. A stator nozzle for a small 120 kW sCO 2 radial inflow turbine has been optimised with a previously developed modularised geometry optimiser based on the Nelder-Mead algorithm. Basic geometry and targets are generated using the preliminary design code TOPGEN. To parametrise the geometry, a parametrised turbine stator mesh generator based on a meridional streamline is developed and used to create a high-quality structure mesh. Then bulk CFD simulations are carried out with the open-source CFD software OpenFOAM. The simulation results are evaluated and a Pareto front is created for optimisation. By applying the developed geometry optimiser, three final different stator geometries have been selected via setting different combinations of weighting factors. They are the optimised Mach number distribution (M σ ), the optimised outlet flow angle distribution (α σ ) and a compromise between the M σ and the α σ cases. The geometries, Mach number contours, normalised total pressure contours and the properties on the outlet boundary patch for these stators are extracted and discussed. These three different optimised stators are as follows.
• The optimised M σ case has the shortest divergent nozzle, returns the best Mach number distribution and performed less loss compared to the baseline case. However, the outlet flow uniformity is not enhanced significantly. • The optimised α σ case has the longest divergent nozzle and returns a more uniform outlet flow compared to the baseline case. However, the Mach number distribution and the efficiency are not enhanced significantly. • The balanced performance stator has a medium length of divergent nozzle compared to the previous two optimised stators. The performances in Mach number distribution and outlet flow angle distribution are also a compromise between the previous two nozzles, which means that both the Mach number distribution and the α are improved, but are not as large as in the optimised M σ and α σ cases.
Both of them, the optimised M σ case, the optimised α σ and the balanced performance stator, are selected as 'good' cases, which can provide good agreement with the targets for a small 120 kW sCO 2 radial inflow turbine. Considering this, the convergent-divergent nozzle and design methodology outlined in this study is the preferred approach for stators. However, there are still some limitations of this study; for example, the Nelder-Mead algorithm is not a global optimisation algorithm, which limits the optimiser to a local minimum. The mesh generation process can only generate structure mesh, which may result in low-quality mesh in some cases, especially when the inclination angle is too large. These limitations suggest a future work using a more advanced optimisation algorithm, such as an artificial intelligence aided algorithm, to optimise the stator geometry. In all, this work will benefit the development and application of the sCO 2 power cycle in the near future.

Appendices Appendix A. Stator geometry parametrisation methodology
To use the optimiser so as automatically to control the generation of the stator blade geometry, the stator blade geometry should be properly parametrised. In the first instance, the optimiser controls the throat radius, width and angle, which sets the position of the floating control points. Next, as discussed Figure A1. Two-dimensional stator nozzle blade parametrisation. Figure A2. Schematic to determine a control point.
in Section 3.2, the transonic stator blade is composed of three different Bézier curves. Thus parametrisation is applied along the lines that control the nozzle shapes. In this appendix, a detailed illustration is shown of how to parametrise the nozzle geometry. The outside part of the divergent section of the nozzle is formed by line 1, which is a Bézier curve with seven control points, as shown in Figure A1. The detailed schematic of line 1 and how to determine a control point is shown in Figure A2. As shown in Figure A1, point A 1&3 is the intersection point between the line XR1-XR2 and the control circle '1' (shown as a green dashed line). Similarly, G 1 is on the trailing edge circle, sharing the co-tangent line with point A 1 , as shown in Figure A3. As shown in Figure A2, points A 1 and B 1 are on the control circle, to control the expansion shape directly downstream of the nozzle throat. The position of B 1 and the length of the arc A1B1 are determined by the angle θ 0 . Similarly, the position of F 1 and the length of the arc G1F1 are determined by the angle θ 1 . Point C 1 is located on the extension of line segment A1B1.
The distance from C 1 to point A 1 is determined by the fraction C1 fA of the length of line A1G1 (shown as the gray dashed line in Figure A2). The term 'fA' denotes the 'fraction along the line'. Hence the position of point C 1 can be determined by the angle θ 0 and the fraction C1 fA . Although the position of point E 1 is defined as similar to point C 1 , the position of E 1 is fixed to reduce the complexity of the optimisation. Point D 1 can move freely along two axes to adjust the curve shape. One axis is along the line A1G1, and another axis is perpendicular to the line A1G1. First, the point D 1 is defined on line A1G1 with the factor D1 fA , which means the length of line A1X1 is D1 fA times the length of line segment A1G1. Then, using point D 1 as the foot of the perpendicular line to find the point D 1 , with the length of line X1 X1 being D1 fV times the length of line A1G1. Line 2, linking A 2 to F 2 , is generated in a similar way. Line 3 forming the outside of the nozzle, has 12 control points.
In engineering applications, even though a sharp corner is possible to make by machining, however its structural strength is poor. Thus, in this project, a finite radius circle is put in the trailing edge section of the turbine stator blade to increase the structural strength. A zoomed-in view of the trailing edge section is shown in Figure A3. As the trailing edge curve affects the wake, these curves can be modified by the optimiser. The position of points G 1 and G 2 can be adjusted with θ 1 and θ 2 .
All of the geometric parameters necessary to define the blade profile (stator throat and trailing edge radii, outlet angle of the nozzle centerline, location of the center of the trailing edge circle, radii of the control circles, etc.) as well as the control points for the Beźier curves can be used as design variables in the optimisation process. Since the converging part of a transonic nozzle only has a minor influence on the flow, the present work is only focused on the optimisation of the diverging part of the blade profile. Fourteen parameters are selected to be adjusted by the optimiser, as shown in Table A1. In Table A1, the subscript 'fA' means 'fraction of the lengths A1G1 or A2G2 along the lines A1G1 or A2G2'; similarly, the subscript 'fV' means 'fraction of the lengths A1G1 or A2G2 vertical to lines A1G1 or A2G2'.
The parametrised stator geometry is defined by a state vector of 14 components, which leads to a 14-dimensional optimisation problem: The modularised geometry optimiser will adjust this vector to find the minimum of the objective function.

Appendix B. Dense gas flow over a 2D backward ramp
At thermodynamic conditions close to the critical point, some dense gases exhibit non-ideal behavior. In this section, to validate the solver's ability to predict non-ideal gas fluid dynamics, a simulation of the non-ideal fluid dynamics of dense gases compared to analytical solutions is presented. It should be noted that this appendix is a concise summary of the study by Qi, Xu, Han, and Jahn (2022). The whole validation process for the solver is published in Qi, Xu, Han, and Jahn (2022). The angle to find point B 3 θ B3C3 The angle to find point C 3 X3 fA The position along the line C3K3 for point X 3 X3 fV The vertical distance to the linee C3K3 for point X 3 R e Exit radius for the calculation domain Cramer (1991) published a detailed study of the non-ideal dynamics of gases, revealing new phenomena including the formation and propagation of expansion shocks, sonic shocks, double sonic shocks and shock splitting, and gaving analytical solutions for these new phenomena. When discussing non-ideal behavior, the fundamental derivative, , is used to determine whether the fluid properties enter the non-ideal gas region (Thompson, 1971). The fundamental derivative is given by

B.1 Theory of non-ideal gas region
where v denotes the specific volume, and the acoustic speed is given by For a perfect gas (p v = R T), it can be shown that Obviously, > 1, if γ > 1 (B4) Thompson (1971) has shown that the region of inverted gas dynamics is defined by M 2 > 1/(1 − ). Here, the Mach number decreases with increasing flow velocity in a steady isentropic expansion. Thus, a region is introduced of non-ideal behavior delimited by J > 0, where J is defined by the relationship between and M: For certain conditions and thermodynamic models, gas can enter this non-ideal region. In this region, classical gas dynamics are inverted, meaning that M decreases through expansion waves as the Prandtl-Meyer function (ν) is inverted. This phenomenon has attracted attention from researchers working on Organic Rankine Cycle (ORC) or sCO 2 cycle design. It is important for a non-ideal gas CFD solver to simulate the fluid dynamics accurately when fluid properties are in this particular region. Durá Galiana et al. (2016) performed a numerical study to show that dense gas properties can enter the non-ideal region through an isentropic expansion from given stagnation conditions. In such an expansion, the Mach number first increases but then decreases, while the fundamental derivative ( ) first decreases to near zero before increasing again to near one. Thus, there is a maximum M during the isentropic expansion, and the non-ideal gas dynamics for 0 < < 1 can be captured. As the expansion process is isentropic and M > 1, the properties along a streamline are described by the analytical solution of an isentropic process.

B.2 Analytical solution
Investigating a similar case and comparing RGDFoam predictions with the analytical solution is a further way to assess the capability of RGDFoam. The fluid considered is MD 4 M, with critical properties of T cr = 653.2 K and P cr = 0.877 MPa. An isentropic expansion routine is designed to determine the analytical solution for the non-ideal expansion process. To ensure that the fluid properties enter the non-ideal region through expansion, the stagnation temperature and pressure are chosen, based on the study by Durá Galiana et al. (2016), to be T 0 = 1.025 T cr and p 0 = 2.0 P cr . A pressure p = 0.001 P cr is applied at the outlet. The analytical solution for the properties along this isentropic expansion process is calculated using properties obtained from the National Institute of Standards and Technology (NIST) real gas database (REFPROP) (Colonna et al., 2008;Lemmon et al., 2013). Using the isentropic assumption, M, T, p and for each point along a streamline through the expansion are evaluated. Figure B1 shows values of p, T and versus local M along the isentropic expansion process. It can be seen that, at the beginning of the expansion, the local M increases as p and T decrease (Figures B1(a) and B1(b)). However, once the expansion enters the non-ideal region, M starts to decrease again. Thus, a maximum local M of 1.962 is reached.

B.3 Simulation case description
For validation of the CFD solver RGDFoam, the fluid dynamics near the local peak in M are of interest. Thus a test case that allows the fluid to expand from the classical region into the non-ideal region is designed. For this, a backward ramp with angle 30 • is selected to expand and accelerate the fluid. Owing to its relative simplicity, the 2D backward ramp is one of the most popular geometries used to evaluate CFD problems. The simulation details are listed in Table B1. Four grids have been chosen for the grid dependency study, with resolutions of 28k, 48k, 75k and 131k. The result is shown in Figure B2. Even though the 48k grid can return a grid-independent result, considering that this case is used to Figure B2. Grid dependency study for the 2D backward ramp.   validate the solver with an analytical simulation, the highest resolution of grid 131k is chosen to carry out the following simulations.

B.4 Results and discussion
Contours of M are plotted in Figure B3 and corresponding properties along two streamlines are shown in Figure B4. Fluid enters the domain with a constant M, 1.8, at conditions corresponding to an isentropic process from the total conditions. Once the fluid reaches the Mach waves originating from the corner, the fluid starts to accelerate. M first increases, while still in the classical region (J < 0), forming an expansion fan. However, once the fluid properties enter the non-ideal region (J > 0), M starts to reduce, signifying a non-ideal process.
To gain a better understanding of this phenomenon, two streamlines are selected from the calculation domain, marked with '1' and '2'. Properties along these streamlines are shown in Figure B4. A close-up view of the M contours along streamline 2 during the expansion process, marked with a rectangular frame in Figure B3(a), is shown in Figure B3(b). Figures B4(a) and B4(b) show the M along the x-position for streamlines 1 and 2, respectively. It is clear from Figure B4(a) that M gradually increases once the fluid starts to turn. For streamline 1, a maximum M is reached at x = 0.32. After this point, the streamlines enter the non-ideal region and the M starts to drop.
Once the streamline is parallel with the second wall, the M is constant again and equal to 1.46. Comparing streamlines 1 and 2, it can be seen that the expansion process is more spaced out, owing to the increased distance from the corner. The maximum and minimum M remain the same, 1.96 and 1.46, respectively.
It can be seen from Figure B4(b) that M increases from point a to point b. An expansion fan is formed as the fluid accelerates and approaches the non-ideal region. Point b marks the position where the maximum M is reached and also the start of the non-ideal flow, as shown in Figure B5. As properties enter the non-ideal region, M starts to decrease.
Considering Figure B4(b), the most rapid reduction in M exists between points c and d. Here, the spatial rate of change of M, ∂M/∂x, is highest for the entire process between points b and e. The region c to d also corresponds to the lowest . This is explained by the local ratio of ν to M ( ν/ M). The values of ν/ M for bc, cd and de are −0.90, −0.69 and −1.46, respectively. Thus, in region cd, the most rapid M decrease per unit turning angle is obtained. For fluids or conditions that result in a lower , this could lead to the formation of a rarefaction shock. Figures B4(c) to B4(f) show the variation of p and T for both analytical solutions and CFD simulations obtained using RGDFoam along the streamlines. It is clear that the pressure and temperature changes against M along both streamline 1 and 2 (bold line) show good agreement with the analytical solution obtained from the MD 4 M expansion routine (dashed line). The implemented RGDFoam solver can accurately recreate the analytical solution. This further demonstrates the ability of RGDFoam to perform predictions for dense gases with operating conditions close to the critical point and in the non-ideal gas region.