Development and application of a modularized geometry optimizer for future supercritical CO2 turbomachinery optimization

ABSTRACT During the operation of supercritical CO (sCO ) turbomachineries, large pressure ratios, supersonic conditions and non-ideal gas fluid dynamic phenomenon may happen, which will decrease the whole cycle efficiencies. Hence non-standard designs for the turbomachineries geometry are needed. Successfully simulations in capturing non-ideal gas fluid dynamics and coupled with optimization algorithm are hard and currently not available in any commercial CFD software. Therefore, the problem of optimizing sCO turbomachinery is difficult to solve, which has become a major obstacle to obtain compact and high-efficiency sCO turbomachineries. In this study, a modularized geometry optimizer is developed to obtain the non-standard geometric designs for sCO turbomachineries. Multiple techniques are applied to this optimizer, include the Nelder–Mead algorithm, Mahalanobis distance and stochastic algorithm. The newly developed optimizer can successfully find the optimum satisfying the objective function under given weighting factors. The computational cost can be reduced through a stochastic algorithm. To validate the optimizer, a convergent–divergent nozzle for air with a target Mach number equal to 2.4 is optimized. Different starting points and combinations of weighting factors are used to create a Pareto front. Adjust the weighting factors for different terms of the objective function leading the optimizer to go to different directions in n-dimensional spaces. Three optimized cases, one is Mach number optimized, one is outlet flow uniformity optimized and the other is compromised case, are picked out and analyzed. The results show that the optimizer can successfully find optimized geometry than the reference case and potentially save computational cost. Due to the modularized characteristics, the components of this optimizer can be replaced with any available techniques, which mean the optimizer can be applied to solve different types of optimization problems.


Introduction
Problems associated with the fast development of human society, such as climate change, environmental pollution, and global warming, have created an international drive toward renewable energy. Among the renewable energy sources, solar energy is the most abundant and it has already been used by human society in various forms since ancient times.
Currently, significant focus is on the development of concentrated solar thermal (CST) energy systems using the supercritical CO 2 (sCO 2 ) Brayton cycles, aiming to a better use for solar energy. Supercritical CO 2 (sCO 2 ) is regarded as an excellent working fluid due to an easier reaching 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 the rapid changes in density reduce the compressor and compressor performance tests to study a 500-kW sCO 2 power generation system for waste heat recovery. Today, the South West Research Institute in collaboration with General Electric are building a test facility to test a 10-MW axial turbine (Kalra et al., 2014).
Much of this axial turbine development is aimed towards future large-scale energy productions ( > 50 MW), where sCO 2 can replace steam in the mid to long term. However, sCO 2 also has a potential for smaller scales (Zhou et al., 2018). For power cycles in the range of 0.1 to 25 MW, using sCO 2 allows a paradigm shift to using efficient radial inflow turbines (RITs). This is the case due to the combined effects of the highly-dense working fluid and comparably low flow rates, resulting in highly power-dense machines.
To acquire a higher efficiency, a good design methodology is required for the sCO 2 RITs. Currently, sCO 2 RITs or turbine components are designed using quasi-onedimensional inviscid or two-dimensional (2D) viscous flow solvers, coupled to the thermodynamics models available at the time (Pasquale et al., 2013;Qi et al., 2017;Zhou et al., 2018). The obtained geometries are based on standard procedures, which are usually standard geometries, i.e. can only apply straight lines or simple curves in blade shape generations. However, risks may happen if only rely on these design procedures. For the RITs working with high pressure ratio, supersonic flow conditions may occur, which will result in shock waves that decrease the efficiency. Or compressors working near the critical points, large fluctuations in densities and heat capacities decrease the stability of operations. These phenomena can hardly be predicted by the quasi-onedimensional inviscid or two-dimensional viscous flow solvers within the traditional design procedures. Thus non-standard geometries are required for the sCO 2 turbomachinery design to minimize the losses due to the above-mentioned effects. A good way to enhance turbines performance is to adjust the three-dimensional (3D) blade geometry. Due to the complexities of the flow and the 3D nature of the passages modifying the shape is a challenging undertaking especially as the optimal shape may be non-intuitive. Hence, it is very hard to get an optimized turbine blade shape.
An alternative approach for the design is the intensive use of computational fluid dynamics (CFD) with appropriate optimization strategies (Abadi et al., 2020;Ghalandari et al., 2019). Usually, people optimized the turbine geometries with adjusting a limited number of variables that cannot really reach to an optimal design. Hence, the research on optimization methods has attracted more attention these years. The aerodynamic shape optimization method can be divided into gradient-based method and non-gradient-based method. The former one relies on the explicit correlations between the objective function and optimized variables. However, for sCO 2 RITs, the correlations are implicit. Thus a non-gradientbased method is more applicable in solving turbomachinery optimization problems, for its nature of no need to explore the relationship between optimization variables and objective function.
As a fact that many attentions were focused on the optimization of sCO 2 power cycle itself, less efforts are paid on the optimization of sCO 2 turbomachineries. Nevertheless, the successful simulations in capturing non-ideal gas fluid dynamics based Riemann problems confronted in large pressure-ratio sCO 2 turbomachineries are very hard, as they need advances methodologies that not available in current commercial CFD software. Hence, it is almost impossible for any available optimizer to automatically manipulate the turbine geometries, update the meshes, operate the non-ideal gas based CFD simulations and evaluate the results simultaneously. Therefore, the problem of truly optimizing sCO 2 turbomachinery has been difficult to solve, which has become a major obstacle to obtain of compact and high-efficiency sCO 2 turbomachinery.
Hence, to overcome these obstacles, in this study, an efficient modularized geometry optimizer is developed, which is based on the coupling of different tools, namely, an optimizer, a geometry generator, a mesh generator, a CFD solver, and a post-processor. The geometry optimizer is designed by modularization methodology, which means different parts of this optimizer can be replaced or changed to any other available techniques. During the searching of optimums, many candidate geometries are analyzed in a fully automated manner. Then the optimizer calls a high fidelity CFD solvers (here in this study we used open-source CFD solver OpenFOAM) to simulation of the given geometries. Based on the results returned by CFD solvers, the optimizer finds the optimal geometry. One more thing to be considered is time consumption during the application of non-gradient-based optimization algorithm. Techniques focused on reducing computational cost should be applied. Thus the objectives of this study are: (1) illustrating the methodology of the modularized geometry optimizer and applied techniques; (2) verification of the developed modularized geometry optimizer; (3) exhibiting an application example in engineering reality.
The structure of this study is listed as follows. Section 2 describes the fundamental methodologies for the Nelder-Mead based optimizer, which includes a description of the Nelder-Mead method and a method for reducing the computational cost. Section 3 presents the validation of the developed optimizer, using an example of optimizing a convergent-divergent air nozzle. Finally, Section 4 gives a summary and future works for this study.

Acronyms cdf
Cumulative distribution function. Greek Symbols The mean of a set of observations. μ The mean or expectation of the distribution. The target value (objective functions) of the optimization progress. φ Variables need to be optimized. σ Standard deviation. ϕ The cost value of rte objective function.

Roman Symbols
Ma Average Mach number [-]. The increment step. t The target value. 4 Radial turbine inlet section. cost Indicates this is a cost value.

Optimization methodology
In this section, the methodology for optimization progress is introduced. This includes a detailed description of the optimization method and a strategy to reduce the computational load. For an optimization problem, we can always define an objective function with a form of f (x 1 , x 2 , . . . , x n ) whose value should be minimized.
To solve this optimization problem, an explicit function solver is applied. However, in many engineering applications, the value of the problem functions we want to minimize may not have explicit relationships to the input parameters. Especially for the turbine design, it is very hard, even impossible to find an explicit relationship between the turbine geometries parameters and turbine performance variables. In that condition, an optimization method that is independent of the explicit control function is needed to minimize the given objective function. Thus this capability requires a modularized strategy to allow us to choose the optimization methodology, which will be discussed in the following sections.

The modularized optimization strategy
The modularized optimization methodology is developed, and the flow chart is shown in Figure 1.As shown in Figure 1, the modularized optimization algorithm consists of four different parts: input/output section, optimizer, function evaluation, and database. The core parts of this modularized optimization methodology are the optimizer and the function evaluation part, shown in the blue and red dashed frames. At the first, the optimizer reads the input files, which includes the control file, initial conditions, ranges, to correctly start the optimization process. Once the optimizer generates a new state vector (will be discussed in Section 2.3), the corresponding objective function is evaluated with a given method (e.g. high fidelity CFD simulation with OpenFOAM or other CFD software). The value of the objective function is then passed to the optimizer for determination, at which part a new state vector may be generated. Once the results satisfied the converging criteria of the optimizer, the output results are provided to indicate the right optimum. During the whole progress, the state vectors and the values of corresponding objective functions are stored in a database, which helps the restart of the optimization progress. Both of the optimizer and the function evaluation modules can be applied in different methods, such as Fletcher and Reeves (Fletcher & Reeves, 1964), Newton, Nelder-Mead, etc. for the optimizer module, or Computational Fluids Dynamics (CFD), Finite Element Method (FEM), etc. for the function evaluation module. The selection of the modules is based on the specific engineering problems. In this study, the future engineering problems are focused on the optimization of sCO 2 RITs, thus, the Nelder-Mead optimization and CFD methods are applied to these two modules. Details about the optimization process will be discussed in the following subsections.

State vector and objective function
For engineering optimization problems, for example, usually, optimize the efficiency through modifying the geometry or the operational parameters. But how to turn the geometry into mathematical parameters, which could be recognized by the optimizer, should be considered. A way to pass geometry and operational information to the optimizer is the parametrization of geometry and then creating state vectors X. The state vector is like this: where V denotes the geometric or operational parameters, for a radial turbine, such as inlet radius r 4 , inlet Mach number M 4 , etc. The state vectors are assigned to the optimizer for evaluation through the simplex transformation method. The state vectors are input variables, then the value of the target of the optimization progress can be denoted as , as In engineering applications, there is typically more than one objective function that needs to be minimized. This results in a multi-objective optimization problem.
In mathematical terms, a multi-objective optimization problem can be formulated as where φ can be replaced to efficiency, power output or temperature, etc. In this study, as the formula for the objective functions is ambiguous, we use a uniformed objective function , defined as where ϕ is the cost value for the objective function, returned from the function evaluation, and W is the weighting factor for different objectives (costs). The sum of the weighted costs of the objective functions is the total objective function , whose value is minimized by the optimizer. One thing to be reminded of is that in multi-objective optimization, there does not always exist a feasible solution that minimizes all objective function values. Hence, the Pareto optimal solutions, i.e. solutions that cannot be improved in any of the objectives without degrading at least one of the other objectives, are used to determine the results of the optimization process. The searching directions of the optimal can be adjusted by setting the weighting factors.

The Nelder-Mead methods
Usually, optimization methods required the local gradient or the first derivative of the objective function to find the optimal solutions, such as the method presented by Fletcher and Reeves (Fletcher & Reeves, 1964), the method presented by Davidon (1991). For turbomachinery or similar problems, it is hard, even impossible to find the local gradient or the first derivative of the objective function, thus derivative-free optimization methods that do not use derivative information in the classical sense to find optimal solutions are required. The Nelder-Mead simplex transform method (Nelder & Mead, 1965) is one of these methods, which is first discussed by J. A. Nelder and R. Mead (Nelder & Mead, 1965) in 1965. The Nelder-Mead algorithm is designed to solve the unconstrained optimization problem by minimizing the value of a given function. It is a simplex-based method that works on a given simplex. This method uses only function values at some point in R n and does not try to extract or calculate an approximate gradient or local derivative at any of the given points. Hence it is one of the general class of direct search methods (Wright, 1996) and is well suited for the turbomachinery blade geometry optimization problems.
A simplex-based method usually starts with a given set of n + 1 points [x 0 , . . . , x n ] ∈ R n . These points are considered as the vertices of the current working simplex S, with the corresponding set of function values at the vertices f j = f (x j ), for j = 0, . . . , n. The method then performs a sequence of transformations of the working simplex S according to the procedure describes in the literature (Singer & Nelder, 2009), aimed at decreasing the function values corresponding to its vertices. At each step, the value of one or more points on the simplex is calculated, together with their corresponding function value. The transformation for the working simplex is determined by comparing the newly calculated value with the other vertices value. This process will be continued till to the working simplex S becomes sufficiently small, or the function values f j are close enough. Thus this algorithm is a 'naturally' time-saving optimization technique.
A suitable initial step allows the optimizer to work more globally. A too-small initial simplex can lead to a local search, consequently the Nelder-Mead method can get more easily stuck in a local minimum. In this project, the initial simplex is determined by the nature of the problem and to make sure a large enough initial simplex is applied to the optimization problem. The Nelder-Mead method needs an N × (N + 1) initial simplex, to make sure that variations for every component are considered, the initial matrix is designed as In this matrix, every single line is a state vector, which consists of N variables, and s denote the increment step. The increments of every variable create a convex hull in the n-dimensional space. The value of increments should be determined according to the given problem.
Once the simplex is introduced, the whole process of the Nelder-Mead method can be suited into the blue frame in Figure 1, and the flow chart for the Nelder-Mead method is shown in Figure 2.To be noted, each optimized variable has a valid range, to eliminate the cases without physical meanings. The optimizer uses a limitation function to enforce this, and the upper and lower limits are read from the input files for each optimization variable. The newly generated state vector is checked to confirm that no variable is outside the range. If the new state vector generated from the transformation of the working simplex is located outside of the range, a high penalty value is assigned to the optimizer to force the optimizer to search in another direction.

Strategy to reduce computational load
Normally, the newly generated state vector is directly passed to the function evaluation module to evaluate the cost, shown as the red frame in Figure 1, which usually consists of a CFD solver. In an optimization context, any gain in computational efficiency for function evaluations is advantageous, especially when using timedemanding CFD simulations to evaluate the objective functions. It can be noticed that the optimizer need many iterations to find an optimum, especially when the simplex arrives at the adjacent zoom of the optimum. Thus it is essential to find a good way to reduce the number of function evaluations. During the operating of the optimizer, an interesting phenomenon is observed. In the earlier stage of optimization progress, the future function evaluations fall outside the convex hull of the current simplex, which leads the working simplex transforming and 'walking' to the optimum. Once the working simplex gets close to the optimum, future function evaluations will fall within the convex hull of the current simplex. Instead of large steps, small steps are then required to converge to the optimum. Usually, the future function evaluation falls within the convex hull of the current simplex and previously visited points. Evaluating these intermediary iterations consumes much computational time. These intermediary function evaluations are almost 'useless' once the optimum is found or nearly found. It is only essential to do an evaluation for obtaining the objective function when the future state is located outside the current working simplex. Hence, the intermediary iterations can be substituted to reduce computational costs. One idea is to replace evaluation with interpolation based on the evaluated data that already exist. In this study, the n-dimensional linear interpolation method is applied to reduce computational cost and enhance the computational efficiency. Figure 3 shows the evaluation method corresponding to the red frame in Figure 1.which includes the interpolation methodology and the high fidelity CFD simulation tool (the green frame). First, the optimizer reads the database and pre-loads the data points, to create a convex hull (which will be discussed in the following section). The newly generated state vector, also considered as a 'point' in n-dimensional space, is first compared to the generated convex hull to determine whether to do interpolation or a new evaluation with the timeconsuming CFD simulation. In order to carry out the n-dimensional linear interpolation, we need to determine if the new state vector is inside or outside the convex hull formed by existing vector points. To avoid the accumulation of errors that may arise from the interpolation of interpolated data, only function evaluations are stored in the database. The interpolation methodology will be discussed in detail in the following sections.
If the state vector needs to be evaluated with CFD simulation, the high fidelity CFD solver will be applied, as shown in the green frame of Figure 3. The high fidelity simulation consists of three parts: pre-processing, CFD simulation, and post-processing. Currently, for this study, the CFD simulation is carried out with OpenFOAM. The values of the objective function include interpolated and evaluated values, are calculated and returned to the Nelder-Mead optimizer. Simultaneously, the state vector, the results from post-processing and the value of objective functions are stored in the database.

Determination of interpolation or evaluation
As shown in Figure 3, there are three 'if' statements to determine whether do interpolation or not. At the starting of the evaluation process, the evaluated points (state vectors, can be treated as a point in n-dimensional space) from the database are ordered by the distance to the new point/state vector. Here we assume these points are only 2D, i.e. only two input parameters. According to the distance to the new point in the R 2 space, two cut off distances R 0 and R 1 are defined, as shown in Figure 3. Different methods are applied, depending on the cut off distance R 0 and R 1 : If the new point N is close enough to a point whose value of objective function already evaluated, use the objective function value of the nearest point f (P n ) as the value of point N.
If the new point N is located between R 0 and R 1 from the nearest point, then an interpolation method is applied. However, an randomize strategy is applied here to guarantee the accuracy of interpolation, will be discussed later.  . Two-dimensional schematic for different methods used in interpolation. (a) N is located within the convex hull and within the distance R 0 from the nearest point. (b) N is located outside of the convex hull but within the distance R 0 from the nearest point. (c) N is located inside of the convex hull and the distance to the nearest point is between R 0 and R 1 . (d) N is located outside of the convex hull and the distance to the nearest point is between R 0 and R 1 . (e) N is located inside of the convex hull but the distance to the nearest point is larger than R 1 and (f) N is located outside of the convex hull and the distance to the nearest point is larger than R 1 . evaluation, as shown in Figure 4, and correspond to Figure 3.
A number of possible scenarios are illustrated in Figure 4. Figure 4(a) shows that if the next geometry point is located inside the convex hull and close enough to a previously evaluated point (point a), then the Nearest interpolation method is applied for point N, i.e. f (N) = f (a). For the second case shown in Figure  4(b), even though the next geometry point lies within the circle around point a with a radius of R 0 , point N lies outside of the convex hull, hence the evaluation method is applied. Figure 4(c) illustrates that point N is located between R 0 and R 1 to the nearest point (point a) and located inside of the convex hull, the linear ndimensional interpolation method is applied. Figure 4(d) illustrates that point N is located between R 0 and R 1 to the nearest point (point a), but outside of the convex hull. Hence, do the evaluation. As shown in Figure 4(e) even though the point N is located inside of the convex hull, the distance between the nearest point is larger than R 1 , thus the evaluation method is applied. Similarly, Figure  4(f) shows that point N is located outside the convex hull and the distance to the nearest point (point c) is larger than R 1 , the evaluation method is applied.

Mahalanobis distance
As discussed above, all the state vector points are arranged by Mahalanobis distance D M . That is because the components of the state vectors have different units. Especially for an optimization problem with multiple design parameters in different units, for example, degree ( • ), meter (m), percent (%), results in a highly nonuniform n-dimensional space. Thus a normalization grid method needs to be applied on the state vectors prior to proximity evaluation. To normalize the grid scale, the Mahalanobis distance method (Mahalanobis, 1936) is employed. The Mahalanobis distance of an observation x = (x 1 , x 2 , x 3 , . . . , x N ) T from a set of observations with mean μ = (μ 1 , μ 2 , μ 3 , . . . μ N ) T and covariance matrix S is defined as Mahalanobis distance can also be defined as a dissimilarity measure between two random vectors x and y of the same distribution with the covariance matrix S, The calculation of Mahalanobis distance (D M ) is performed with the spatial.distance.mahalano bis, package from open-source scientific tools for Python . For the initial step, the covariance matrix is formed from the initial simplex. Subsequently, the covariance matrix is formed from all available geometries. Then according to Equation (6), the Mahalanobis Distance is calculated.

Randomization strategy for acceleration
However, once the state vector points are within the range R 0 to R 1 , purely applying interpolation may cause some loss of details, even causing failing to find the real optimum. One way to mitigate this is to randomly select a limited number of state vector points for interpolation. The rest cases are still evaluated with high fidelity CFD simulations. That method provides accurate data for interpolations, guarantees the accuracy of the final optimum.
The accuracy for the interpolation decreases with the distance from the data points used for interpolation. Thus to optimize the randomization interpolation progress, the possibility of interpolation should be increased with decreasing distance to the nearest point. For N within the range of R 0 to R 1 to the nearest point, the evaluation or interpolation is determined randomly selected. A threshold value T r is set, a random number N r is generated, then • if for N r ≥ T r , an evaluation is performed; • else for N r < T r , an interpolation is performed.
In order to adjust the possibility to do interpolation, the threshold value T r should be linked to the distance. A good way is to determine the threshold T r by the cumulative distribution function (cdf ) of a normal distribution.
In probability theory and statistics, the cdf of a random variable X, evaluated at x, is the probability that X will take a value less or equal to x. In the case of a continuous distribution, it gives the area under the probability density function from minus infinity to x. The cumulative distribution function of a normal distribution is where μ shows the value where the cdf is 0.5, the mean of the corresponding standard deviation function. If we set μ = 1 2 (R 0 + R 1 ) that guarantees around 50% state vector points will be calculated through interpolation methods, while the rest 50% will be evaluated through high fidelity CFD simulations. σ is the standard deviation of the cdf. If we set σ = 1 4 (R 1 − R 0 ), i.e. (R 1 − μ) = (μ − R 0 ) = 2σ , the cdf is at 2.5% and 97.5% at R 1 and R 0 .
In the following example, when setting R 0 = 1 × 10 −7 , R 1 = 1 × 10 −2 , the corresponding cdf with 1 2 (R 1 − R 0 ) = 2σ is plotted in Figure 5.In this example, the value of the cdf is increasing from 2.5% to 97.5% as the Mahalanobis distance approaches 0.01. During the determination section, for example, if the D M of the next state vector point N to the nearest point is 0.004, the value of cdf is 0.3452. This means an average 34.5% of cases at this distance will be evaluated. This method gives a continuous transmission from a high possibility of doing interpolations to a high possibility of doing a new evaluation with CFD simulations as D M increases.

The method to create a minimum convex hull in n-dimensional space
In order to carry out the interpolation, it is essential to create a convex hull in the n-dimensional design space. If using all of the data to do the linear interpolation in n-dimensional space, the interpolation time will increase dramatically. A good way to reduce the interpolation time is to only use the closet points. First, the minimum convex hull covering the new geometry point is identified, and then interpolation is performed with the vertices of the minimum convex hull.
The convex hull for the geometry point cloud is found by the Quickhull algorithm (Barber et al., 1996;Harris et al., 2020). First, the convex hull of the whole data cloud is obtained. Next, it is determined, if N is inside the convex hull or not. Next step, the farthest point (determined by Mahalanobis distance D M ) is removed from the point cloud, and the new convex hull is generated. Then again it is assessed if point N is inside the reduced convex hull. This procedure repeats, until the optimizer finds the minimum convex hull that covers the point N in ndimensional space or until the point cloud has reduced to n + 1 points in n-dimensional space.

Interpolation method
There are many different types of interpolation methods, such as linear interpolation, cubic interpolation, Kriging interpolation, etc. This modularized optimizer can also choose different interpolation methods. In this study, linear interpolation is applied to carry out the interpolation in n-dimensional space. This method is from the Python class LinearNDInterpolator of the Scipy package .
After the development of a method to reduce the computational cost, optimization progress is tested with and without the interpolation methods. The results show that the acceleration ability is determined by the time of a single CFD simulation. Once a single CFD simulation use less or equal time to a single interpolation, the acceleration effect is not obvious. In that case, only doing CFD simulation is a wise choice to guarantee the accuracy of the optimization. However, once the optimizer applies heavy CFD simulation that cost much longer time than a single interpolation, the acceleration effect is quite significant.

Validation of the modularized geometry optimizer
This section presents an example to show the capability of the geometry optimizer. What is more important, this example also validates the operation of the geometry optimizer.

Case description
The verification and validation process needs to show that the optimizer can find a geometry that is better than the baseline design, which exhibits conditions of optimality. In this section, a convergent-divergent nozzle with an exit Mach number of 2.4, taken from Ref. Anderson (1990) is used as a reference solution to be studied and optimized. The reference nozzle shape is shown in Figure 6.
The nozzle shown in Figure 6 is a convergent-diver gent air nozzle designed by the method of characteristics. The target outlet flow is a uniform flow with a Mach number of 2.4. For the reference case, the area ratio is 2.33, which is within 3% of the value A e A * = 2.403 obtained for isentropic 1-D Mach number relations. This small error is induced by the graphical construction of the nozzle. This nozzle shape is calculated analytically, which may not return uniform flow with a Mach number of 2.4 in reality. Hence the goal for the CFD optimization in this example is to use the nozzle designed by the method of character and to generate an improved geometry. Improved geometry means that the nozzle has a more uniform flow and that the exit Mach number is closed to the target of 2.4.

Methodology
In this section, the methodology of the optimization of a convergent-divergent air nozzle is presented. Figure 7(a) shows the lines and blocking strategy used for the mesh generation for the nozzle.

Mesh geometry parametrization and blocking strategy
It can be seen from the figure that the inlet section or the subsonic section of the nozzle consists of two blocks, BLK0 and BLK1. The nozzle shape of the inlet section is controlled by line a 0 a 2 , which consists of a straight line a 0 a 1 and a Bézier curve a 1 a 2 . The control points of line a 1 a 2 are c 0 , c 1 and c 2 . The throat is defined by the points a 2 and b 2 (b 2 is located at origin).
As the problem for the transonic nozzle is hyperbolic, which means the downstream will not affect the upstream flow field of the nozzle. To reduce the complexity of the  problem, the geometry of the subsonic section of the nozzle is set to constant. The nozzle throat width is set to a fixed value to keep a constant mass flow rate and to reduce the complexity of the optimization problems.
The supersonic section of the nozzle is one single block, BLK2. The nozzle shape is defined by another Bézier curve, line a 2 a 3 . This curve is defined by six control points (d 0 to d 5 ), to be noted that the number of control points are user defined, here the number 6 is only for illustration. Initially, the six control points are uniformly located along the line a 2 a 3 . Positions of these control points and a3 form the geometry vector, where x and y indicate the x-and y-coordinates of the corresponded points. However, practically, it would be better to write the x-and y-coordinates in a relative form to point a3: Similar, · · · · · · dn y = a2 y + l n (a3 y − a2 y ).
This method allows us having a better imagination for where these points are. Obviously, the control points near the throat will have stronger influence on the downstream flow field than others. Based on the mathematical nature of the Bèzier curve that the first and last control points determine the curve entry and outlet slopes. Hence, to make sure the smoothness along the nozzle throat wall, the slope at the left-end of a2a3 is set to zero with adding a control point near the a2 and set the y-coordinate equal to a2 y . Then the geometry vector can write in another format that: For the current nozzle, the control point away from the nozzle section will have less impact on the downstream flow filed. Thus, to reduce the complexity of the simulation problems, we only use five control points to initialize the nozzle geometry. The examples of initial simplexes are presented in Appendix section. Point a3 determines the nozzle length and nozzle exit area ratio, which is given by The optimizer adjusts the nozzle geometry by modifying the geometry vectors during the optimization progress. For the reference nozzle, the a2a3 is graphically fitted from Figure 6, and the coordinates of key points are presented in Table 1 to help reproducing the simulations. For the baseline case of the optimization problems, the initial simplexes are also presented in Appendix.
The mesh of the convergent-divergent nozzle is created by the geometry package of the in-house code Eilmer4 (Jacobs & Gollan, 2017;Jahn et al., 2018), and the foamMesh utility to convert the mesh to OpenFOAM format. A 2D structured mesh with constant clustering towards the wall is generated, which is shown in Figure 7(b). The 2D problems for the convergent-divergent nozzle are symmetric, thus only half part of the nozzle is simulated.

Simulation set-ups
The inviscid CFD simulation is done with the opensource library OpenFOAM (Foundation, The Open-FOAM, n.d.). The transonic/supersonic, compressible gas transient solver sonicFoam (Foundation, The Open-FOAM, n.d.) is used to solve the RANS equation. Total pressure and total temperature boundary conditions are applied to the inlet boundary, shown in Figure 7(b). The waveTransmissive boundary condition is applied to the nozzle outlet boundary for pressure field to allow the passing of shock waves. The pressure ratio is set to a larger value which guarantees that choke conditions occur at the nozzle throat. A slip wall boundary condition and a symmetry boundary condition are applied to the nozzle wall and center line, respectively. The turbulent models are switched off and the inviscid condition is applied to give a better comparison with the analytical solution. Detailed simulation setups are presented in Table 1.
Three grids with the cell numbers of 513, 2033, 8463 have been used to carry out a grid dependency study. The result of the grid dependency study is shown in Figure 8, which is the value of Mach number from the throat to the nozzle exit along the center line.It can be seen from the figure that the results from the medium size mesh are in good agreement with the result from the fine mesh. Thus the medium size mesh is used in the following simulations.

Objective functions and postprocessing
During the optimization, the post-processor carries out postprocessing procedure and returns the values for cost evaluation.
The optimizer is used to reduce the cost for the objective functions, thus robust and meaningful objective functions are needed. As a nozzle with uniform Mach 2.4 outlet flow should be obtained, the objective function should be focused on two terms, the exit flow Mach number and the exit flow uniformity. To achieve this we select an objective function of the form: Here ϕ total denotes the total value of the objective function, ϕ denotes the cost of every term, and W denotes the weighting factor. The optimizer explores the design space to reduce the total cost. The individual cost for Mach number and flow direction should consider both the target value and the uniformity. This can be achieved by calculating the standard deviation for a given value, for example, the Mach number, through the following equation: where i corresponds to different cell faces across the nozzle exit. If replacing Ma by Ma t , the equation is considering the standard deviation from the target Mach number. As the area is different for faces, the weighting with respect cell area is required. Thus the cost for Mach number is calculated as Similarly, the cost for outlet flow angle is used to considering the uniformity of the outlet flow, which is written  as For an ideal nozzle, the value of Ma cost and α cost is zero. Hence the objective function becomes A good way to carry out optimization progress is to try different combinations of simplex and weighting factors to allow the optimizer to scan the design space. In this example, the starting vectors for the initial simplex are selected randomly. The nozzle length is selected between 6 and 12 times the nozzle throat width, which covers the analytical solution, whose value is 7.9 times of the nozzle throat (Anderson, 1990). Then different combinations of weighting factors are applied to the optimizer to obtain different optimal solutions.

Results and discussion
Different initial simplexes with different combinations of weighting factors result in a total number of 4196 CFD runs. The Appendix section gives four examples for initial simplexes and combinations of weighting factors. With these initial simplexes and weighting factors, the optimizer can do a search in n-dimensional space. Figure 9 shows the Mach cost value versus α cost value The dot corresponded to the example nozzle from the literature (Anderson, 1990) is shown in this figure in red, located at the bottom side. This nozzle has good uniformity of the exit flow, however, the Mach number cost value is not optimal. That is because the area ratio for the analytical solution is about 2.3 rather than 2.403. This figure confirms that the optimizer successfully finds nozzle geometries that are better than the reference nozzle, with more uniform flow and the exit flow Mach number close to the target.
It can be seen in Figure 9 that there are some local 'shining' regions, with a lot of geometries cumulated around these regions, like the adjacent regions of (0.02, 5.7), (0.01, 2.2), (0.04, 4.0), (0.15, 0.5), etc. These regions indicate the local optimal in the case spaces related to different nozzle lengths. That's because the Nelder-Mead optimization algorithm is not a global algorithm, hence we need to apply multiple initial simplexes and weighting factors. For the current case, instead of allowing the optimizer to adjust the nozzle length, different fixed lengths are set to carry out the optimization. Different initial simplexes and weighting factors will affect the locations and the distribution of the dots. With multiple tries, a lot of CFD cases are generated, then form the current figure. However, no matter how we set the initial simplexes or the weighting factors, there is a line the dots cannot lie beyond, which is the Pareto front, formed by the cumulated CFD runs at the left bottom side of the figure. That is because, the ideal nozzle with the target Mach number and a 100% parallel exit flow will have an infinite length, which is practically unreality. Thus for a nozzle with finite length, there is a trade-off between the Mach number and the exit flow uniformity. The Pareto front shows this balance. Figure 10 shows a zoomed in view of Figure 9 for the left bottom part.As shown in Figure 10, three points are selected and marked as Case A, Case B and Case C. Case B has the lowest value of Mach cost, Case C has the lowest value of α cost. Case A is a balance of Mach cost and α cost, which is selected from the Pareto front, and has the minimum distance to the origin (0, 0) of all simulated  geometries. The details for CASE A to CASE C, including points coordinates, simulation setups are presented in Table 1.
The nozzle shape of these three cases and the analytical nozzle geometry are plotted and compared in Figure 11. Case A has the longest nozzle length. Case B, which has the lowest Mach cost, shows a divergent exit at the outlet of the nozzle. It is always best to have a uniform Mach number 2.4 flow with a continuous divergent part of the nozzle. Case C has the most uniform outlet flow compared to the other cases. Figures 12(a-d) show the Mach number contours for different cases. It is obvious that for the analytical nozzle from the literature (Anderson, 1990), the outlet Mach number does not reach 2.4 yet, as shown in Figure 12(a).
For Case A, Mach number 2.4 contour almost occupies the whole outlet section and good uniformity of the exit flow exist. For Case B, which is shown in Figure 12(c), even though the nozzle exit flow is closer to the target number, the uniformity is low, as shown by the flow streamline. Case C is the nozzle with the best exit flow uniformity, as shown in Figure 12(d), however, the exit Mach number is away from the target Mach number 2.4.
To have a better comparison of the outlet flow for these nozzles, Figures 13(a,b) show profiles of exit flow Mach number and flow angles distribution. Case B has the best distribution of Mach number, and the area average Mach number is almost equal to the target Mach number 2.4. The reference nozzle from the literature (Anderson, 1990) does return an exit flow whose Mach number is less than 2.3. This confirms that the optimizer can truly find a geometry that returns an exit flow with almost zero Mach deviation cost, as shown by Case B. Figure 13(b) illustrates the distribution of flow angles α for the four nozzles. It can be seen from the figure that Case B has the worst distribution of flow angles, especially near the wall, here the exit flow angle is almost 8.0 • . To have a better understanding of the flow angle distribution, vectors have been plotted on the data dots. It can be seen that Case C has the best uniformity for the outlet flow. These results confirm that the optimizer finds various optima, depending on the objective function.
One extra item to be added into this context is the optimization convergence history. Figure 14 shows the convergence history for the optimization process and different state vector components for Case B.It can be noticed from Figure 14(a) that large fluctuation in total cost happens in the first few steps. That is due to the initializing of the initial simplex, which is large. After about 50 steps, the simplex starts to shrink on to approach the local minimum. Then after about 150 CFD runs, the state vector no longer changes, which implies that the local optimal is found.
To be reminded that the optimizer uses the acceleration technique which is described in the previous section. The red line exhibits the running history for optimization without acceleration technique. After 200 runs, the optimization without acceleration technique starts to shrink on a target value. However, this optimization process is still far away from the local minima. Hence, this figure provides evidence that the randomization algorithm reduced the computational time significantly. Figure 14(b) shows the changing of each geometry parameters corresponding to the CFD runs. It confirms that after about 50 steps the geometry parameters are tending to be steady, and after about 150 CFD runs, the geometry parameters no longer change, which implies the optimum is found.

Conclusion
In this study, a modularized geometry optimizer is developed. Multiple techniques are applied to this optimizer, which includes the Nelder-Mead simplex algorithm, Mahalanobis distance algorithm, randomization acceleration algorithm and linear interpolation methods. The application of these techniques can potentially reduce the computational cost of the optimization iterations. The new developed modularized geometry optimizer has the following highlights: (1) Can successfully find the optimum satisfying the objective function (2) Through adjusting the weighting factors, different optimal can be attained (3) A randomization algorithm, supported by the interpolation method is implemented to reduce computational cost (4) Optimizer is verified through the optimization of a convergent-divergent air nozzle (5) Developed with modularization strategy that allows future continued development and application to other problems.
To test the optimizer, a convergent-divergent nozzle for air with a target Mach number equal to 2.4, is optimized. Different starting points and a combination of weighting factors are used to create a Pareto front. The results show that the optimizer can successfully find optimal geometry over a baseline case. Three different objectives have been chosen and analyzed. One is the lowest Mach number cost, and one is the lowest exit flow angle cost, the final one is a compromise between the two objectives, and has the lowest total cost. This shows that adjusting the weighting factor for different terms of the objective function causes the optimizer to go in different directions to return different results. The optimizer is suitable to be applied to other fluid dynamics problems.
However, there are still some limitations for this piece of work, such as that the Nelder-Mead algorithm is not a global one, thus we need to have multiple tries to avoid stuck in local optimum; lack a criteria for choosing initial simplexes; need a more stable method to calculate the convex-hull, etc. In future work, the modularized optimizer will be used to optimize the radial inflow turbine blades for sCO 2 utilization.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This section provides the initial simplexes to gain the CFD results.