Numerical optimisation of a micro-wave rotor turbine using a quasi-two-dimensional CFD model and a hybrid algorithm

Abstract

Wave rotors are unsteady flow machines that exchange energy through pressure waves. This has the potential for enhancing efficiency over a wide spectrum of applications, ranging from gas turbine topping cycles to pressure-gain combustors. This paper introduces an aerodynamic shape optimisation of a power generating non-axial micro-wave rotor turbine and seeks to enhance the shaft power output while preserving the wave rotor’s capacity to function as a pressure-exchanging device. The optimisation considers six parameters including rotor shape profile, wall thickness, and number of channels and is done using a hybrid genetic algorithm that couples an evolutionary algorithm with a surrogate model. The underlying numerical model is based on a transient, reduced-order, quasi-two dimensional computational fluid dynamics model at a fixed operating condition. The numerical results from the quasi-two-dimensional optimisation indicate that the best candidate design increases shaft power by a factor of 1.78 and imply a trade-off relationship between torque generation and pressure exchange capabilities. Further evaluation of the optimised design using three-dimensional computational fluid dynamics simulations confirms the increase in power output at the cost of increased entropy production. It is further disclosed that increased incidence losses during the initial opening of the channel to the high-pressure inlet duct compromise the shock strength of the primary shock wave and account for the decrease in pressure ratio. Finally, the numerical trends are validated using experimental data.

1 Introduction

1.1 Operating principles and characteristics

A wave rotor dynamically exchanges energy between gas streams of high enthalpy and gas flows of low enthalpy through means of shock and expansion waves. Opposed to turbomachinery, where the energy is transferred aerodynamically, wave rotors pressurise through wave action. Shock wave compression features high efficiencies and can take place over small volume making pressure-exchange machinery particularly attractive for small-scale power applications. The geometry of a wave rotor does not feature any complex surfaces and comprises simple channel shapes. Due to the presence of both the hot expanded gas stream and the cold air stream within the same machinery, the rotor is inherently self-cooled and the material temperature remains below the peak cycle temperature.

Fig. 1

Principal wave rotor components of a four-port device consisting of wave rotor arranged circumstantially and flanked to each side by inlet and outlet stator endplates

Several layouts of wave rotors have been introduced ranging from radial [1, 2] and non-axial [3,4,5,6] to axial [7] and from throughflow [8], where the fluid predominantly flows in one direction only, to reverse flow models [9]. In addition, depending on the type of application the number of ports per cycle can vary between two-port [10], three-port [11], four-port [12], five-port [13], and nine-port [14] devices. The present paper is dedicated to a throughflow four-port configuration with non-axial channel shape. In such a design, the wave rotor consists of a cylindrical drum with cambered channels arranged around the circumference. The rotor is located in between two stator end plates, which house discrete manifold openings. An example of such a wave rotor with curved channels is shown in Fig. 1. The basic principle of wave rotors revolves around energy exchange through shock and rarefaction waves. In the example design shown in Fig. 1, it is the aim of using wave action to transfer energy from the high-pressure gas (HPG) to pressurise the low-pressure charge induced via the low-pressure air (LPA) duct. When rotating about its axis, the rotor channels are periodically and almost instantaneously subjected to the pressure and temperature in the manifolds. This sparks the generation of shock and expansion waves. In addition, cambered wall profiles in combination with angled port manifolds introduce momentum change and enable the wave rotor to be applied as both a pressure-exchange device and a power turbine.

Fig. 2

a Schematic of wave rotor turbine arrangement within a gas turbine as one possible application for wave rotor turbines. The concept features a four-port throughflow wave rotor with cambered channels for shaft power extraction and is connected to a combustor through the high-pressure zone (HPG-HPA). b Unwrapped view of the wave rotor in the ${{\theta }}$–z plane outlining distribution and location of shock and expansion waves

The main arrangement of manifolds for a four-port throughflow device is shown in the schematic of Fig. 2a. Heat addition may take place using a steady-flow combustor as done in gas turbines or an internal combustion engine. The ducts comprise a high-pressure inlet and outlet as well as a low-pressure inlet and outlet. To examine the behaviour within the channels during operation, it is helpful to follow a single channel as it rotates around the circumference and unfold its path on to a two-dimensional ${{\theta }}$–z plane, as done in Fig. 2b. This essentially allows to determine the distribution and location of shock and expansion fans as well as the distribution of hot and cold gas streams. The arrangement of ports needs to be carefully timed with the arrival of shock and expansion waves in order to maximise the strength of shock waves, limit attenuation, and reduce backflow of gas from the channels into the inlet ducts.

Initially, air enters the wave rotor at low pressure and temperature through the low-pressure air (LPA) at station 1. Compression of this low-pressure charge air is achieved via three shock waves (S1, S2, and S3). The two strongest compression waves (S1 and S2) are generated at the bottom of Fig. 2b when the channels are exposed to the outlet from the combustor, namely the high-pressure gas (HPG, station 3). Prior to this, the gas within the channels are at pressures close to ambient. However, as soon as the channel is exposed to the high-pressure conditions prevalent within the HPG port (station 3), a primary right-travelling shock wave (S1) is triggered. As soon as S1 reaches the right-hand-side end of the channel, the high-pressure air (HPA) outlet is opened and the shock wave is reflected, resulting in a secondary left-running shock wave (S2). Primary and secondary shock waves essentially exchange energy from the high-pressure gas coming from the combustor to the low-pressure charge air before it exits the wave rotor at an elevated pressure level through the HPA duct (station 2).

Closing the HPG port as soon as the secondary shock wave impinges on the inlet side effectively stops the flow and inhibits spillage of gas back from the channel into the duct. In addition, port closure initiates an expansion or rarefaction wave (R1). This, in combination with an additional expansion fan (R2) that is generated upon exposing the channels to the LPG duct, gradually expands the hot exhaust gases towards the LPG port. Compared with a traditional gas turbine arrangement, rarefaction waves R1 and R2 effectively represent the function of a turbine. Finally, a weaker shock wave (S3) is the result of closure of the low-pressure gas (LPG, station 4) port.

The ratio of inlet mass flow rates between HPG and LPA, termed loop flow ratio ($\lambda ={\dot{m}}_{\mathrm {HPG}}/{\dot{m}}_{\mathrm {LPA}}$), mainly governs scavenging of the ports of hot exhaust gases and thus the cold flow penetration length $L_{{\mathrm {p,c}}}$. In general, $\lambda$ will be greater than unity and $L_{{\mathrm {p,c}}}<L$, which results in a certain amount of exhaust gas recirculation (EGR) being looped back to the combustor inlet. Furthermore, if the hot gas penetration length $L_{{\mathrm {p,h}}}$ does not suffice to push the fresh charge air towards the HPA outlet, fresh air exhaustion (FAE) takes place and some of the compressed air gets expelled unused.

1.2 Literature survey

Wave action devices have been subject to interest from a variety of research institutions. Earlier investigations started by the Brown Boveri Company in the 1940s sought to apply wave rotors as a topping device in gas turbines aiming at increasing system efficiency [15]. The concept has been further developed by Rolls Royce [16,17,18], ONERA [19, 20], and NASA [21,22,23]. In addition to gas turbines, wave rotors have successfully been applied as pressure wave superchargers to both light- and heavy-duty internal combustion engines [16, 24,25,26]. Recent research studies have taken up the concept and applied wave rotor technology to refrigeration cycles [27,28,29] and pressure-gain combustion [10, 30,31,32,33]. Wave rotors characterised by a non-axial channel shape to extract shaft power while operating as a pressure-exchange device have first been used in a gas turbine arrangement by Pearson [34]. The device achieved a power output of approximately 26 kW with a thermal efficiency of 10%. The wave rotor unit comprised helically shaped passages and angled port ducts to maximise momentum transfer. Recently, Michigan State University and Warsaw University of Technology have developed a wave-disk engine concept housing a radial wave rotor unit with curved channels yielding thermal efficiencies of up to 10% [35].

Early wave rotor studies resorted to graphical methods, such as the method of characteristics, to design and analyse the compressible wave action within wave rotors. With the advent of computational fluid dynamics (CFD) and rising computational power, elaborate one- and two-dimensional codes have been developed in order to examine wave rotor designs with greater accuracy, such as those done by Paxson at NASA Glenn Research Center for axial wave rotors [36,37,38,39,40]. The model includes source terms to account for leakage and wall friction and account for finite passage opening effects through modification of the boundary conditions. The code was further used to examine the dynamic behaviour [41] and area variation effects [42] and was further extended to account for channel curvature through a passage averaged approach [40]. Further one-dimensional codes were developed at the Mathematical Science Northwest Inc. [43], ONERA [19, 20], and the Warsaw University of Technology in both one and two dimensions [44, 45].

Compared to the number of one- and two-dimensional studies, there are relatively few three-dimensional simulation studies available in the open literature. This is based on the fact that they are by far more resource intensive. Piechna et al. [46] compared the results of two-dimensional with three-dimensional results and found—besides a highly skewed contact surface—that centripetal and Coriolis accelerations cause a distortion of moving shock waves.

There has been no attempt in the open literature to optimise the wave rotor shape in order to maximise performance. The requirement to examine a sufficiently large number of design candidates is exacerbated by the need for transient CFD simulations for the evolution of moving shock and expansion waves. A quasi-steady approach using a moving frame of reference (MRF) is not feasible. This paper, for the first time to the author’s knowledge, thus seeks to introduce a reduced-order CFD model used to optimise the channel shape of an axial throughflow wave rotor turbine. To reduce the computational cost, a hybridised genetic algorithm (GA) with a systematic design space exploration and exploitation strategy was used. The study will demonstrate that the reduced-order model in combination with the hybrid optimisation method can efficiently and swiftly deliver an optimised wave rotor shape and reproduce the trends witnessed in experiments and higher-order three-dimensional CFD. The structure of the paper is thus as follows. Section 2.1 introduces the baseline wave rotor model and compares its design with other existing wave rotor models. Section 2.2 presents the scope of the optimisation study and the main principles of the hybrid algorithm, followed by the numerical set-up and domain discretisation for the quasi-two-dimensional model along with results from the mesh and numerics sensitivity study. In Sect. 2.2, the numerical set-up of the three-dimensional CFD model is given, while Sect. 2.4 gives a brief overview of the experimental layout. The results of the optimisation study are discussed in Sect. 3.1, which are further discussed by a loss and gas dynamic analysis using 3D CFD. Finally, the trends witnessed in the numerical models are compared to experimental data.

2 Methodology

2.1 Baseline wave rotor

The wave rotor turbine forming the basis for this study and used as a baseline model for the optimisation has been extensively experimentally tested [47] and used for Q1D validation [6]. The prototype features symmetrically cambered wall profiles, following the shape of a parabola. A computer-aided design (CAD) image of the baseline model is shown in Fig. 3 illustrating the curved wall profile with further details on the geometric dimensions of the rotor given in Table 1.

Fig. 3

Baseline geometry with symmetrically cambered channels

Table 1 Geometric dimensions of baseline wave rotor turbine

In order to compare different wave rotor designs, one can employ the relations developed by Wilson and Fronek [48] and Nagashima et al. [49], who defined non-dimensional parameters for finite passage opening T, viscosity F, and leakage flow G. These are given as

$$\begin{aligned} \begin{aligned} T&= \frac{W_{{\mathrm {ch}}} a}{U_{{{\theta }}} L}\\ F&= \frac{\sqrt{\nu L/a}}{D_{\mathrm {h}}}\\ G&= \frac{2\delta }{h} \end{aligned} \end{aligned}$$

(1)

where $W_{{\mathrm {ch}}}$ denotes channel width, a the speed of sound, $U_{{{\theta }}}$ tangential velocity, and L the channel length. In addition, $\nu$, $D_{{\mathrm {h}}}$, $\delta$, and h are the kinematic viscosity, hydraulic diameter, axial leakage gap, and channel height, respectively. This signifies that finite passage opening time is directly proportional to passage width and speed of sound, while being inversely proportional to the tangential velocity and channel length. Thus, the smaller the channel width or the greater the channel length and rotational speed, the smaller the finite opening effects. Frictional effects can be expected to be greater for channels of greater length and smaller hydraulic diameter. Finally, leakage is directly related to the leakage gap, as well as channel height.

A comparison of these non-dimensional parameters for a variety of wave rotors from the open literature is provided in Table 2 along with geometric dimensions and rotational speed. Although most of the listed proponents are designed to be used purely for pressure exchange and were tested at relatively low temperatures, it becomes apparent that the baseline model in this study features similar characteristics in terms of friction and finite opening timing effects. Leakage in the present model is variable and comparable to ABB’s Comprex and the University of Tokyo’s design. It is noteworthy that the figure holds only for cold conditions. In reality, however, thermal expansion will alter the figure. Compared with previous designs, the presented baseline model is considerably smaller in diameter, length, and channel height. Nonetheless, in terms of mean tip speed the rotor is within comparable range.

Table 2 Comparison of previous wave rotors, both actually realised ones that were tested and conceptual ones (University of Tokyo/ONERA) with the wave rotor model used in this study (adapted from [49])

2.2 Optimisation campaign

2.2.1 Wave rotor optimisation

The optimisation seeks to achieve increased torque generation through modifications to the wall camber profile. This approach maintains the inlet and outlet stator port solution, as well as channel/port height and the overall length of the rotor. In addition, the number of design parameters is reduced as well with benefits regarding computational effort, albeit at the cost of limiting its capacity for exploitation. With respect to experimental validation, it further ensures that the same stator end plates can be employed and merely the rotor itself needs to be replaced rather than the entire wave rotor assembly for back-to-back comparisons.

To facilitate automatised geometry updates, the wall camber profiles are approximated via Bézier curves with five control points giving full control over the contour shape, as shown in Fig. 4. Additional parameters are the number of channels, which, in combination with channel width, dictate the wall thickness $t_{\mathrm {w}}$. Therefore, the optimisation deals with six design parameters in total. Lower and upper bounds to each of the parameters then establish the design space and constrain the optimisation (Table 3). The first control point is located at the leading edge and remains fixed. The subsequent control points are fixed with respect to the horizontal location, but are allowed to move in azimuthal direction.

Fig. 4

Bézier curve parameterisation of the wave rotor wall camber

Table 3 Lower and upper bounds for design variables for the optimisation. It is noteworthy that for ease of implementation into the optimisation routine, the azimuthal bounds for control points along the channel camberline were defined in terms of length rather than angles

The aim is to maximise the time-averaged power output ${\overline{P}}$ while maintaining the pressure-exchange capabilities of the wave rotor. In this study, the energy transfer capacity and thus efficiency are estimated through the achieved pressure ratio. To this end, the objective function contains one main function aimed at maximising power and three additional constraints for designs that do not meet the criteria specified for pressure ratio, channel width, and the number of channels. It is formulated as

$$\begin{aligned}&\mathrm{OF}(x_i) = {\left\{ \begin{array}{ll} 10 \bigg | \frac{P_{{\mathrm {t}}} - {\overline{P}}}{P_{{\mathrm {t}}}} \bigg |, &{} {\mathrm {if}}\ {\overline{\varPi }}_{{\mathrm {TT}}} \ge 0.97 {\overline{\varPi }}_{{\mathrm {TT,BL}}} \\ 10-{\overline{\varPi }}_{{\mathrm {TT}}} &{} {\mathrm {if}}\ {\overline{\varPi }}_{{\mathrm {TT}}}< 0.97 {\overline{\varPi }}_{{\mathrm {TT,BL}}}\\ 20-1000 |W_{{\mathrm {ch}}}-0.0023| &{} {\mathrm {if}}\ W_{{\mathrm {ch}}} < 0.0023\ \text {and}\\ &{} N_{{\mathrm {ch}}}\%n_{\mathrm {cyc}}\ne 0 \\ 20 &{} {\mathrm {if}}\ N_{{\mathrm {ch}}}\%{{ n _{\mathrm {cyc}}}}=0 \\ \end{array}\right. }\nonumber \\ \end{aligned}$$

(2)

where $P_{\mathrm {t}}$ and ${\overline{P}}$ represent a predefined target power and the estimated average power output for a given candidate, respectively, while ${\overline{\varPi }}_{{\mathrm {TT}}}$ denotes time-averaged pressure ratio, $W_{{\mathrm {ch}}}$ channel width, $N_{{\mathrm {ch}}}$ the number of channels, and $n_{{\mathrm {cyc}}}$ the number of cycles per revolution. The main objective function is defined as a weighted relative error from the estimated power output to a target power output $P_{\mathrm {t}}$. The weighting factor was selected as 10. The optimiser then attempts to minimise the error between the actual power output and the target power output and in the process favours designs with a higher power output.

The first additional constraints for energy transfer capacity and manufacturability are defined. Energy transfer through the shock waves may be compromised by mixing losses and incidence losses as the flow enters the channel from the ports. Therefore, one can expect a trade-off relationship between the ability to generate more power, while functioning as a pressure-exchange device at the same time. Thus, in order to allow the optimiser to seek those designs with maximum power output and enable as many designs as possible to be considered without being discarded, a threshold is defined that permits a 3% decrease in compression pressure ratio compared to the baseline model. This is reflected by the factor 0.97. The fitness value of each design candidate surpassing this limit is judged based solely on the predicted power output. The higher the power output, the better the objective function value of the design candidate. If, however, a particular design fails to achieve a compression ratio higher than the defined threshold, a high penalty is assigned. This decreases in a linear fashion the higher the pressure ratio and the closer it is to the predefined threshold.

In the same fashion as the baseline model, design candidates are required to be manufacturable through machining from solid. This translates into a minimum channel width of 2.3 mm, which is reflected in the third constraint within the objective function. If the channel width of a candidate design is smaller than the 2.3-mm threshold, then a linear penalty function is assigned with a weighting factor of 1000 attributed to the difference between the candidate channel width and the minimum channel width. Finally, as power generation in wave rotors is generally pulsatile [47], the final optimisation constraint function within the objective function aims at generating power in a as smooth as possible manner. This is done by imposing an offset between subsequent cycles and thus avoiding large torque amplitudes. In detail, this is achieved by assigning a large penalty to cases where the total number of channels is dividable by the number of cycles without any remainder. In such a scenario, the channel arrangement in each cycle would exhibit the same exact relative position with respect to the ports and would thus exacerbate torque amplitudes.

Fig. 5

Flow chart for transient wave rotor optimisation coupling an optimisation routine with CFD

The steps taken during the optimisation are summarised in the flow chart displayed in Fig. 5. The optimisation starts with the user parameterising the geometry and selecting the corresponding bounds. These are then normalised in the MATLAB routine to a hypercube. In addition, the objective function is set along with the maximum number of expensive CFD computations. It is noteworthy that the optimisations can be continued from a given data set should it transpire that the initial computational budget was too optimistic. As soon as this is done, the actual optimisation is started following the framework outlined in Fig. 6. At each iteration, the updated geometry is fed into the CFD solver and updated and a new mesh created. Then, the CFD simulation is run based on the baseline flow field in order to minimise the number of iterations to reach a limit cycle, which is based on convergence of mass as well as pressures and temperatures in the ports.

Finally, an averaging period is appended to obtain time-averaged data for predicted power ${\overline{P}}$, given as

$$\begin{aligned} {\overline{P}}= \frac{1}{t_{\mathrm {end}}{-}t_0} \int _{t_0}^{t_{\mathrm {end}}} \frac{\tau (t) 2\pi N}{60} \mathrm{d}t \end{aligned}$$

(3)

where $\tau (t)$ represents instantaneous torque, t time, and N the rotational speed. Time-averaged pressure ratio ${\overline{\varPi }}_{{\mathrm {TT}}}$ is evaluated by

$$\begin{aligned} {\overline{\varPi }}_{{\mathrm {TT}}}= \frac{1}{t_{\mathrm {end}}{-}t_0} \int _{t_0}^{t_{\mathrm {end}}}\varPi _{{\mathrm {TT}}}(t) \mathrm{d}t \end{aligned}$$

(4)

where $\varPi _{{\mathrm {TT}}}(t)$ denotes the instantaneous total-to-total pressure ratio. In both instances, instantaneous properties are integrated across the averaging period defined by starting time $t_0$ and end simulation end time $t_{\mathrm {end}}$. Subsequent division through this time difference results in the required time-averaged properties, which are fed into the optimisation routine.

2.2.2 Surrogate modelling

In general, optimisation problems can be mathematically expressed as

$$\begin{aligned} \begin{array}{ll} {\mathrm {Minimise}}&{}\quad f(x_i)\quad {\mathrm {with}}\ i=1,2,\ldots ,n;\\ {\mathrm {Subject\ to:}} &{}\quad x_i^{{\mathrm {lb}}} \le x_i \le x_i^{{\mathrm {ub}}},\\ &{}\quad g(x_i) \le 0 \end{array} \end{aligned}$$

(5)

where the objective function $f(x_i)$ usually remains unknown. In addition, the corresponding design parameters $x_i$ constitute the design space and are subject to lower and upper bounds, while the problem dimension is dictated by n. Furthermore, the objective function can be subject to additional constraints $g(x_i)$. In case of the wave rotor optimisation, evaluating $f(x_i)$ requires elaborate unsteady RANS simulations. This signifies considerable computational effort, in particular when considering several hundreds of these expensive simulations are necessary to cover the design space.

As a result, it would be advantageous to carry out the optimisation on a function that is considerably cheaper (i.e., quicker) to execute. This can be done through replacing the original costly objective function $f(x_i)$ with an auxiliary problem resulting in

$$\begin{aligned} \begin{array}{ll} {\mathrm {Minimise}}&{}\quad {\hat{f}}(x_i)\quad {\mathrm {with}}\ i =1,2,\ldots ,n;\\ {\mathrm {Subject\ to:}}&{}\quad x_i^{{\mathrm {lb}}} \le x_i \le x_i^{{\mathrm {ub}}},\\ &{}\quad {\hat{g}}(x_i) \le 0 \end{array} \end{aligned}$$

(6)

where the surrogate models ${\hat{f}}(x_i)$ and ${\hat{g}}(x_i)$ replace the expensive original objective and constraint functions $f(x_i)$ and $g(x_i)$. The surrogate models are an approximation of the actual objective function, but are significantly cheaper to solve. Ideally, one would seek to be able to extract all the necessary information from the surrogate model without the necessity to run the expensive black-box function $f(x_i)$. However, before this can be done, it is necessary to train the surrogate model. In this paper, the surrogate model is constructed using a Kriging or Gaussian process regression method implemented in the DACE [50] library. The predicted output by the Kriging model,

$$\begin{aligned} \begin{aligned} {\hat{f}}(x_i) = \mu (x) + Z(x) = \sum _{i=1}^{k} y_i(x)\beta _i + Z(x) \end{aligned} \end{aligned}$$

(7)

comprises a deterministic term $\mu (x)$, realised with polynomial function, and a stochastic term Z(x), represented by a Gaussian random function. The stochastic term Z(x) represents the error between the true model and the surrogate model and is approximated through a Gaussian correlation function with zero mean and constant variance. While DACE offers a number of different choices for regression and correlation, the two selected ones appeared to perform best when applied to the test cases introduced in Appendices 1 and 2. The basis function and the corresponding regression coefficients are given by $y_i(x)$ and $\beta _i$, respectively. The regression model is represented by ${\hat{f}}(x_i)$, which can follow zeroth-, first-, or second-order polynomial forms. In the course of this study, the regression term is modelled through a second-order polynomial which follows the form ([51]):

$$\begin{aligned} \begin{aligned} \mu (x) = \beta _0 + \sum _{i=1}^{k} x_i \beta _i + \sum _{i=1}^{k}\sum _{j=1}^{k} \beta _{ij} x_i x_j \end{aligned} \end{aligned}$$

(8)

where the regression coefficients $\beta$ can be determined using least-square methods minimising the error between the regression model and the actual function.

Fig. 6

Hybrid algorithm flow chart

2.2.3 Hybrid genetic algorithm

The first step in the optimisation routine is to generate a set of training data. In most cases,, for surrogate optimisation this is done via design of experiments, such as Latin-hypercube sampling [52, 53]. For the purpose of this study, we propose to run the expensive objective function using evolutionary algorithms, such as a genetic algorithm (GA) or particle swarm optimisation (PSO), to produce the initial training set. After this is accomplished, the routine proceeds to build the Kriging surrogate model on which the optimisation routine is continued. This represents an approximation of the costly original problem. The estimated optimum predicted by the surrogate model optimisation is subsequently compared against its objective function value by running the expensive black-box function. The information from that evaluation is then used to update the surrogate model. This measure is necessary as for most multidimensional engineering problems, the initial surrogate model will not be detailed enough to infer an optimal combination of design parameters.

Following the suggestion of Regis and Shoemaker [54], an iterative procedure is applied that alternates between local exploitation and exploration. The proposed method seeks to determine the minimum of the auxiliary problem while balancing both exploratory and local search. The former signifies that the algorithm seeks to search regions that have not yet been explored well, while the latter aims at further exploring and exploiting reasonably well-explored regions in an effort to isolate a local optimum.

In detail, this signifies that one needs to determine the point within the design space $\varOmega$ with the maximum possible distance to any previously evaluated points as

$$\begin{aligned} \begin{aligned} {\tilde{x}} = \max \limits _{x \in \varOmega } \min \limits _{1\le i\le n} |x - x_i | \end{aligned} \end{aligned}$$

(9)

This point is termed as the maximin point. Hence, by defining a search pattern $\beta$ vector one can set the subsequent evaluation point to be at least $\beta {\tilde{x}}$ from all previously evaluated design points. If the subsequent evaluation point is in close proximity to previously evaluated points, the search is more geared towards local exploitation. Conversely, a point that is far from all previously evaluated points facilitates global exploration of the design space. This is implemented in an alternating manner and included as an additional constraint for solving the auxiliary optimisation problem. Figure 6 summarises the main steps of the hybrid optimisation routine, starting from the initial training set provided by a GA, which is run for a prespecified number of generations, which is dependent on the dimensionality of the optimisation problem (usually in the range of 2–5 generations). Based on the initial data set, a Kriging surrogate model is created before iterations commence solving the auxiliary optimisation problem periodically switching between exploration and exploitation. At the end of each iteration, the expensive problem is run and the surrogate model subsequently updated. As it is not possible to ascertain whether a minimum is local or global, the optimisation is run until the computational budget is exhausted and the maximum number of expensive evaluations reached.

Fig. 7

a Example of the Ackley multimodal two-dimensional test function. b Convergence and evolution of the surrogate model employing the hybrid GA on the Ackley test function at three different instances throughout the optimisation

The principle of the optimisation algorithm is exemplified in Fig. 7 when applied to the Ackley multimodal test function with two parameters. The approximation of original function by the surrogate model at three different instances during the optimisation run is given. The first one after 51 iterations displays the surrogate mode just after completion of the initial GA. It further shows that the next data point to be evaluated is located at the extremities of the search space as a consequence of a large $\beta$ value of 0.95. After approximately 157 iterations and despite the large number of local minima, the routine is swiftly able to determine the location of the global minimum near (0, 0). Further continuation of the search algorithm results in further exploitation of the region in the close vicinity of the global minimum as well as exploration of behaviour further away from it.

In order to further facilitate the search of the design space, a normalisation of the parameter ranges to a hypercube has been implemented. It is expressed as

$$\begin{aligned} \begin{aligned} x_{{\mathrm {norm}}} = \frac{x_i - {\mathrm {min}}(x)}{{\mathrm {max}}(x) - {\mathrm {min}}(x)} \end{aligned} \end{aligned}$$

(10)

In addition, objective functions characterised by a large spread between minimum and maximum function values are normalised using the paired log transformation introduced by Regis and Shoemaker [55]. It is expressed as

$$\begin{aligned} \text {plog}(x)={\left\{ \begin{array}{ll} \text {log}(1+x) &{} \quad \text {if }x \ge 0\\ -\text {log}(1-x) &{} \quad \text {if }x < 0 \end{array}\right. } \end{aligned}$$

(11)

Applying this transformation to the data set alleviates the discrepancy between extreme values, while retaining the location of maxima and minima. This is illustrated in the one-dimensional Rastrigin function in Fig. 8. In order to avoid an already flat to flatten out further, the paired log transformation is only applied if the maximum detected value after the initial GA data set is three times higher than the median of all previously evaluated values.

Fig. 8

Effect of paired log transformation on the one-dimensional Rastrigin function

Table 4 Summary of boundary conditions used for the mesh sensitivity study

2.3 Numerical campaign

2.3.1 Quasi-two-dimensional CFD model

In order to evaluate the performance of different wave rotor designs, transient simulations were conducted using the commercial CFD solver Ansys Fluent R19.1. The governing equations for continuity, momentum, and energy were solved using the unsteady RANS equations with turbulence closure provided through Menter’s k-$\omega$ SST model with enhanced wall treatment. Thus, the flow is assumed to be fully turbulent, which is justified based on the wave rotor channel and port dimensions of approximately 2–5 mm and 15 mm, respectively, as well as the mass flow rate under investigation (${\dot{m}}_{{\mathrm {HPG}}}=32\,\hbox {g/s}$). This results in Reynolds numbers in excess of $1.5\times 10^6$.

Wave rotor simulations are inherently transient despite steady-state conditions prevailing in the ports. Throughout the optimisation study, the boundary conditions set in the inflow and outflow ports are kept constant. The high-pressure inlet port (HPG) is designated as a mass flow inlet with a prescribed inlet mass flow rate of 32 g/s and 773 K inlet temperature. The static pressure in the corresponding high-pressure outlet port is controlled via a user-defined function (UDF) in Fluent such that the mass flow error to the high-pressure port inlet mass flow rate (i.e., 32 g/s) is minimised. As far as the low-pressure region is concerned, the inlet mass flow rate in the LPA duct is again controlled via a UDF, resulting in a loop flow ratio ${\dot{m}}_{{\mathrm {HPG}}}/{\dot{m}}_{{\mathrm {LPA}}}$ of 1.7 entering at 300 K. The low-pressure outlet port (LPG) features a static backpressure according to ambient conditions (Table 4).

The computational domain encompasses the rotor passages as well as the inlet and outlet stator regions, which are connected via a sliding mesh interface, as shown in Fig. 9. Since the approach aims at reducing the computational cost, the mid-plane was extracted at 60 mm between hub and shroud rendering the domain two-dimensional in the z–$\theta$ plane. The stator regions contain the constant axial leakage gap of 0.3 mm on inlet and outlet sides, respectively, and feature three cycles consisting of a high-pressure inlet/outlet and a low-pressure inlet/outlet. Thus, the model does, however, not contain the leakage cavities and also cannot account for windage losses of the rotor within its casing.

Fig. 9

Computational domain for the quasi-two-dimensional optimisations containing stationary inlet and outlet stator domains and a rotating region containing the rotor passages

Fig. 10

Domain discretisation using hexahedral elements with a single cell in radial direction

Table 5 Mesh sensitivity study and corresponding resolution for each sub-domain

Domain discretisation was done using hexahedral elements in Ansys Mesher. As the general flow direction follows the positive z-axis, this approach ensures reduced deviation of face normals and the flow vector and thus limits the amount of numerical diffusion. Furthermore, the discretisation employs a single cell in radial direction with each side being modelled as a symmetry surface. An example of the discretised domain and mesh elements is shown in Fig. 10. This renders the approach quasi-two-dimensional. Prior to the optimisation study, a mesh convergence study was conducted. In total, four different mesh resolutions were evaluated ranging from approximately 58,000 to 350,000 cells, as presented in Table 5.

Hybrid initialisation was used for all simulations and the simulations run for 2500 time steps with a time step of $10^{-6}\,\hbox {s}$ with a maximum number of 15 iterations and a CFL number of 25. The choice of time step corresponds to a rotation of $0.2^\circ$ in between time steps. This was followed by an averaging period for pressure ratio and predicted output power for period of 300 iterations. In addition, all simulations were run employing the implicit, coupled density-based solver in double precision mode. The operating medium was assumed to be air behaving as an ideal gas with dynamic viscosity–temperature dependency approximated by Sutherland’s law. With respect to the numerical set-up, second-order schemes were used for inviscid fluxes and turbulent kinetic energy and dissipation rate, while an implicit first-order formulation was used for temporal discretisation.

Fig. 11

Normalised results for averaged predicted power and pressure ratio as a function of cell count. Data from the finest grid solution were used for normalisation

Table 6 Results of studies on the effect of numerical schemes for inviscid fluxes, turbulence, and transient formulation

The results of the grid sensitivity study are shown in Fig. 11. The graphs depict the averaged predicted power and pressure ratio against the number of cells within the domain. The data were normalised by the solution of the finest grid and reveal that the error between the coarsest and the finest solutions for power and pressure ratio is 5.5% and 2%, respectively. For power, there is virtually no change between medium and fine mesh solutions, while for pressure ratio coarse, medium, and fine mesh solutions are well within 1%. Looking at the elapsed run time in Table 5 shows that there is a significant cost–benefit associated with running on a medium or coarse mesh resolution. In the light of the large number of CFD evaluations necessary in the optimisation, it was decided to proceed the optimisation study with the coarse mesh and accept a slightly larger numerical uncertainty with regard to pressure ratio and power at the benefit of rendering the simulations substantially cheaper.

Fig. 12

Domain discretisation in three dimensions using hexahedral elements for inlet stator, outlet stator, and rotor channels

The next step is to look into the sensitivity of the result with respect to the numerical schemes, the time step size, and the number of inner iterations. The results of this study are summarised in Table 6. At first, the turbulence scheme was set to first-order upwind, which resulted in a 4.1% lower predicted power and a 1.5% higher pressure ratio. Computational effort remains more or less constant. Switching to first-order upwind scheme leads to a larger deviation in terms of predicted power, while pressure ratio remains within 2.5%. The advantage with regard to CPU cost is as expected due to expedited convergence. Reverting to second-order flow and turbulence schemes and switching transient formulation to second-order results in a 3.7% predicted power reduction and 2% increase in pressure ratio, while slightly increasing CPU costs. Finally, a simulation has been carried out with all numerical schemes set to second-order upwind with the exception of spatial discretisation being approximated through third-order MUSCL scheme. This results in moderate changes to the predicted averaged power and pressure ratio, but a 10% increase in computational effort.

The subsequent study deals with the effect of time step size on the numerical accuracy. All studies were run with 15 as a maximum number of iterations per time step. Continuously increasing time step size from approximately $0.12^\circ$ per time step up to $0.2^\circ$ per time step does not seem to impair the solution accuracy and enhances computational efficiency. Approaching $0.48^\circ$ or higher appears to have a more distinct effect on the recorded solution and exhibits a penalty on convergence rate.

Finally, the effect of the number of maximum iterations per time step was examined. This was varied between 20 and 7. Throughout this study, it showed that increasing this parameter to the maximum did not result in a significant change in the convergence rate and residual level at the end of each inner iteration. The residuals continued to decrease approximately two orders of magnitude within seven iterations, thus achieving the set relative convergence target of 0.01.

Based on the findings of the quasi-two-dimensional study, it was decided to proceed with the coarse mesh solution, second-order spatial discretisation scheme for convection and turbulence terms, while using first-order temporal discretisation. All simulations are carried out with a time step of $1\times 10^{-6}$ s and a maximum number of seven inner iterations. Using this set-up, each design point in the optimisation requires approximately 1.8 hours to perform geometry and mesh update and run the simulation on an Intel Xeon E5 3.40 GhZ processor of 12 physical cores and 24 threads.

2.3.2 Three-dimensional CFD model

The quasi-two-dimensional model represents an abstraction and simplification from the actual computational domain aiming at facilitating the computational effort while retaining major characteristics and flow features. In order to gain fidelity in the optimisation results and enhance the comparison between the best candidate design from the pool of design points with the baseline design, further in-depth 3D-CFD simulations are carried out. The objective of this is to include further details in the geometry that were omitted in the quasi-two-dimensional study to reduce computational expense. This includes elements such as inlet guide vanes shown in Fig. 12. The domain itself, however, is similarly structured with two stationary domains housing inlet and outlet as well as a rotating domain representing the rotor. Channel-to-channel interaction is modelled through an axial leakage gap, while leakage cavities and paths around the rotor are omitted due to the complexity and expenditure of such an in-depth transient simulation. Similar to the quasi-two-dimensional approach, the domain is meshed using Ansys ICEM CFD employing hexahedral elements for stationary and rotating domains. The mesh statistics and quality metrics are given in Table 7. The total number of cells used in the simulation was approximately 7.2 million cells. The rotor domain was discretised with 100 cells in longitudinal direction and $16\times 30$ cells for the channel cross section. Overall, the minimum mesh quality throughout the domain is 0.34, minimum face angle approximately $20^\circ$, and maximum skewness 0.8.

Table 7 Mesh statistics of the 3D CFD simulations

The density-based solver of Ansys Fluent R19.1 was used in double-precision mode to solve the unsteady RANS equations with an additional species equation being solved allowing to track the gas entering through the high-pressure inlet and the low-pressure inlet separately and be able to infer the amount of FAE and EGR. Convective fluxes were modelled using Roe’s flux-difference splitting approach with a third-order MUSCL discretisation for the flow equations, while second-order upwind schemes were used for turbulence equations. Temporal discretisation was done using a second-order implicit formulation with a time step size of $1\times 10^{-6}\,\hbox {s}$ and a maximum of ten iterations per time step. All boundary conditions applied are in line with the quasi-two-dimensional study. The working medium was defined as air behaving like a calorically imperfect gas, where both heat capacity at constant pressure $c_p$ and dynamic viscosity $\mu$ are defined as functions of temperature. In case of the former, temperature dependence is defined through a seventh-order polynomial yielding

$$\begin{aligned} \begin{aligned} c_p(T)&= 999.36 - 0.064T + 2.93\times 10^{-4}T^2 + 1.58\times 10^{-7}T^3 \\&\quad - 4.54\times 10^{-10}T^4 + 2.84\times 10^{-10}T^5 - 7.51\times 10^{-10}T^6 \\&\quad + 7.40\times 10^{-21}T^7 \end{aligned}\nonumber \\ \end{aligned}$$

(12)

where T is in K and $c_p(T)$ is in J/(kg$\cdot$K). Dynamic viscosity was modelled following Sutherland’s law.

Along with average pressure ratio and predicted power output, the effects of the optimisation on compression and expansion efficiencies shall be examined. These are determined through considering hot and cold gas streams separately and taking EGR and FAE into account, as presented by Chan et al. [56]. In addition, internal heat transfer is also taken into account as done by Wilson et al. [57]. Therefore, compression and expansion efficiencies are defined as the ratio of actual work to the isentropic work while distinguishing between the compressed, expanded, EGR, and FAE streams. They are given as

$$\begin{aligned} \begin{aligned} \eta _{{\mathrm {comp}}}&= \frac{{\dot{m}}_{{\mathrm {LPA}}} c_{p,{\mathrm {c}}} \big (T_{{\mathrm {HPA}}}^{{\mathrm {is}}} - T_{{\mathrm {LPA}}} \big ) + {\dot{m}}_{{\mathrm {EGR}}} c_{p,{\mathrm {h}}} \big (T_{{\mathrm {HPA,h}}}^{{\mathrm {is}}} - T_{{\mathrm {LPG}}}\big )}{ {\dot{m}}_{{\mathrm {LPA}}} c_{p,{\mathrm {c}}} \big (T_{{\mathrm {HPA,c}}} - T_{{\mathrm {LPA}}} - \Delta T_{{\mathrm {c}}} \big ) + {\dot{m}}_{{\mathrm {EGR}}} c_{p,{\mathrm {h}}} \big (T_{{\mathrm {HPA,h}}} - T_{{\mathrm {LPG}}} - \Delta T_{{\mathrm {c}}} \big ) }\\ \eta _{\mathrm {x}}&= \frac{{\dot{m}}_{{\mathrm {LPG}}} c_{p,{\mathrm {h}}} \big (T_{{\mathrm {HPG}}} - T_{{\mathrm {LPG}}} - \Delta T_{{\mathrm {x}}} \big )+ c_{p,{\mathrm {c}}} {\dot{m}}_{{\mathrm {FAE}}} \big ( T_{{\mathrm {HPA}}} - T_{{\mathrm {LPG,c}}} - \Delta T_{{\mathrm {x}}} \big ) }{{\dot{m}}_{{\mathrm {LPG}}} c_{p,{\mathrm {h}}} \big (T_{{\mathrm {HPG}}} - T_{{\mathrm {LPG}}}^{{\mathrm {is}}} \big ) + c_{p,{\mathrm {c}}} {\dot{m}}_{{\mathrm {FAE}}}\big (T_{{\mathrm {HPA}}} - T_{{\mathrm {LPG,c}}}^{{\mathrm {is}}} \big )}\\ \end{aligned}\nonumber \\ \end{aligned}$$

(13)

where ${\dot{m}}$ denotes the mass flow rates and T the temperature corresponding to a specific port. The additional temperature terms $\Delta T_{\mathrm {c}}$ and $\Delta T_{\mathrm {x}}$ account for the heat transfer between hot and cold gas streams within the wave rotor. These are generally not known and therefore assumed to be equal. In order to be able to solve this equation, compression and expansion efficiencies are assumed to be equal $\eta _{{\mathrm {cx}}}=\eta _{\mathrm {comp}}=\eta _{\mathrm {x}}$.

Local evaluation of irreversibilities and thus potential penalties in efficiency are analysed on the basis of local entropy generation rate. Entropy production is large in regions of large temperature gradients, as is the case in mixing regions of cold and hot gas streams as well as along the hot/cold interface surface. In addition, losses occur in areas of increased turbulent dissipation as well as velocity gradients, such as in shear layers and the boundary layer. According to Iandol and Sciubbia [58] and Herwig et al. [59], the total local entropy generation rate per unit volume ${\dot{S}}_{{\mathrm {p,}}V}$ reads

$$\begin{aligned} \begin{aligned} {\dot{S}}_{{\mathrm {p}},V} = {\dot{S}}_{{\mathrm {p},{\bar{\mathrm {D}}}}} + {\dot{S}}_{{\mathrm {p,D'}}} + {\dot{S}}_{{\mathrm {p},{\bar{\mathrm {C}}}}} + {\dot{S}}_{{\mathrm {p,C'}}} \end{aligned} \end{aligned}$$

(14)

and comprises entropy production terms representing direct viscous dissipation from the averaged flow field (${\dot{S}}_{{\mathrm {p},{\bar{\mathrm {D}}}}}$), turbulent dissipation (${\dot{S}}_{{\mathrm {p,D'}}}$), heat transfer through the average temperature field (${\dot{S}}_{{\mathrm {p},{\bar{\mathrm {C}}}}}$), and heat transfer due to fluctuating temperature gradients (${\dot{S}}_{{\mathrm {p,C'}}}$). In detail, these terms can be defined in Cartesian coordinates as

$$\begin{aligned} {\dot{S}}_{{\mathrm {p},{\bar{\mathrm {D}}}}}= & {} \frac{\mu }{{\overline{T}}} \bigg [2 \bigg [\bigg (\frac{\partial {\overline{u}}}{\partial x} \bigg )^2 + \bigg (\frac{\partial {\overline{v}}}{\partial y} \bigg )^2 + \bigg (\frac{\partial {\overline{w}}}{\partial z} \bigg )^2 \bigg ]\nonumber \\&+ \bigg (\frac{\partial {\overline{u}}}{\partial y} + \frac{\partial {\overline{v}}}{\partial x} \bigg )^2 \nonumber \\&+ \bigg (\frac{\partial {\overline{u}}}{\partial z} + \frac{\partial {\overline{w}}}{\partial x} \bigg )^2 + \bigg (\frac{\partial {\overline{v}}}{\partial z} + \frac{\partial {\overline{w}}}{\partial y} \bigg )^2\bigg ] \nonumber \\ {\dot{S}}_{{\mathrm {p,D'}}}= & {} \rho \frac{\varepsilon }{{\overline{T}}}\nonumber \\ {\dot{S}}_{{\mathrm {p},{\bar{\mathrm {C}}}}}= & {} \frac{\lambda }{{\overline{T}}^2} \bigg [ \bigg ( \frac{\partial {\overline{T}}}{\partial x} \bigg )^2 + \bigg ( \frac{\partial {\overline{T}}}{\partial y} \bigg )^2 + \bigg ( \frac{\partial {\overline{T}}}{\partial z} \bigg )^2 \bigg ]\nonumber \\ {\dot{S}}_{{\mathrm {p,C'}}}= & {} \frac{\mu _{{\mathrm {turb}}} c_{p}}{\mathrm{Pr}_{{\mathrm {turb}}} {\overline{T}}^2} \bigg [ \bigg ( \frac{\partial {\overline{T}}}{\partial x} \bigg )^2 + \bigg ( \frac{\partial {\overline{T}}}{\partial y} \bigg )^2 + \bigg ( \frac{\partial {\overline{T}}}{\partial z} \bigg )^2 \bigg ] \end{aligned}$$

(15)

where $\mu$ denotes dynamic viscosity, T the static temperature, u, v, and w the velocity components in the respective coordinate directions, $\varepsilon$ the turbulent dissipation rate, $c_{p}$ the specific heat at constant pressure, and $\mathrm{Pr}_{{\mathrm {turb}}}$ the turbulent Prandtl number assumed to be 0.9. The equations for entropy generation have been implemented in Ansys Fluent via a custom field function.

Integrating local entropy generation rate ${\dot{S}}_{{\mathrm {p,}}V}$ across a volume yields the total entropy generation rate, which is as follows

$$\begin{aligned} \begin{aligned} {\dot{S}}_{{\mathrm {p}}}&= \int _V {\dot{S}}_{{\mathrm {p}},V} \mathrm{d}V \end{aligned} \end{aligned}$$

(16)

which can be further integrated over a time interval to give the total entropy generation, defined as

$$\begin{aligned} \begin{aligned} S_{\mathrm {p}}&= \int _{t_0}^{t_1} {\dot{S}}_{{\mathrm {p}}} \mathrm{d}t \end{aligned} \end{aligned}$$

(17)

Fig. 13

a The baseline wave rotor with symmetrically shaped channels, b its optimised version. c Schematic showing the experimental layout and sensors used throughout experimental testing

2.4 Experimental set-up

After machining the optimised rotor, both the baseline and optimised rotor were tested experimentally on the University of Bath gas stand. Both rotor designs can be seen inFig. 13a, b. The experimental set-up used for validation follows an open-loop configuration, as shown in Fig. 13c. A detailed explanation of the set-up is beyond the scope of this paper, but follows the same layout as described in [47].

In general, the set-up features capabilities to vary dynamometer excitation load, loop flow ratio $\lambda$, mass flow rates, and inlet pressures for both the high- and low-pressure inflow ducts. In particular, pressurised air is directed to the wave rotor and can be heated up to target inlet temperature through the use of two 44-kW electrical heaters. Loop flow ratio between HPG and LPA is controlled via two pneumatically actuated variable orifice control valves, while backpressure in the HPA duct is also modulated through an additional gate valve on the outlet side. Throughout testing, pressures and temperatures at the wave rotor inlet and outlet are recorded through pressure transducers and k-type thermocouples. Finally, wave rotor speed is controlled via a water-cooled eddy-current dynamometer.

3 Results and discussion

3.1 Wave rotor optimisation and validation

The optimisation routine for the hybrid-GA considered 342 different designs, while MATSuMoTo dealt with 450 different design candidates. Plotting predicted power output versus compression pressure ratio for both optimisation routines indicates a trade-off relationship between power and pressure ratio with a Pareto front, as illustrated in Fig. 14. The results imply that there is a balance between the amount of energy that can be extracted through the change in momentum and the amount of energy being exchanged through the shock and rarefaction waves. The results shown in the graph were normalised by the baseline power output and pressure ratio with the baseline design marked by a red square. The 97% pressure ratio threshold introduced in (2) is illustrated by the dashed line. The best candidate design across both optimisation campaigns, outlined as the filled green circle, is located at this threshold and indicates an increase in power output of 1.78 at the cost of a 3% decrease in pressure ratio. Comparing the two optimisation techniques implies that MATSuMoTo features more of an emphasis on local search, while the results of the hybrid method are more scattered throughout the design space as the algorithm constantly iterates through the predefined search pattern.

Fig. 14

Results from the optimisation results for a the hybrid-GA method and b the MATSuMoTo surrogate model

Fig. 15

a Convergence rate of objective function and normalised power output of hybrid method and MATSuMoTo. b Direct comparison of camberline distributions and Bézier control points for the baseline and optimised geometry

Fig. 16

Numerically determined instantaneous power generation of a single channel as it rotates around the circumference using 3D URANS. The top graph shows the power distribution over multiple cycles showing the pulsatile nature of power generation. The detailed view below compares power distribution of baseline geometry (solid, black line) with the optimised model (dash-dotted, red line) along with cycle-averaged values

However, as can be seen in the convergence plot shown in Fig. 15a both methods swiftly decrease the objective function and find candidate designs with a significant increase in the predicted power output. In this instance, the hybrid method appears to converge sooner after fewer than 100 function evaluations, while MATSuMoTo requires approximately 350 iterations to determine a design with a similar objective function. However, it shall be noted that the final difference in the objective function and thus the estimated power benefit is relatively small. In addition, the optimisation settings (i.e., GA and LHS parameters as well as search patterns) may have the potential to be further optimised for the given optimisation problem. A comparison of the wall contour camberlines is shown in Fig. 15b. As opposed to the baseline design, the optimised design no longer features a symmetric shape. In addition, camber increased with a higher deflection in the centre and more negative channel angles towards $z/L=1$. On top of the geometric modifications, the number of channels has been reduced to 43, resulting in greater channel width.

3.2 Further examination of 3D CFD results

In order to gain further insights and fidelity in the results of the quasi-two-dimensional model, 3D URANS simulations on both baseline and the optimised rotor have been performed. Tracking a single channel as it rotates around the circumference and logging instantaneous power results in the distribution shown in Fig. 16. The top graph shows instantaneous power across multiple cycles (i.e., approximately three rotations). The detailed view below extracts power distribution for a single cycle and compares results for baseline geometry (solid, black line) with the optimised model (dash-dotted, red line) along with cycle-averaged values.

Fig. 17

a Instantaneous entropy generation rate distribution as a function of azimuth angle. b Direct comparison of compression and expansion efficiency $\eta _{\mathrm{cx}}$ and integrated entropy generation within a single channel between baseline and power optimised geometry. The data were normalised by the corresponding baseline values

All data shown were normalised by the average power output of the baseline design. Negative values denote shaft power extraction, while positive values represent power supplied to the rotor. First of all, it becomes apparent that the baseline camberline distribution is far from optimal with respect to extracting work. As expected, the power generation distribution emphasises the pulsatile manner, in which torque is generated in the wave rotor. The main contribution of torque generation stems from the opening of the HPG port to the channels. During approximately two-thirds of the cycle, the channel is idle and does not contribute to power generation. In addition, there are extended periods with positive power, reducing the overall power output.

The optimised model, on the other hand, exhibits increased loading and thus an greater power output. Overall, the increase outlined by 3D CFD is approximately 1.85 and thus confirms the findings from the quasi-2D model used in the optimisation study. Thus, the quasi-2D model is able to yield representative results at a fraction of the computational cost of a full 3D simulation. The reason for the small discrepancy between quasi-2D and 3D simulation results most likely resides in the greater model accuracy, such as the inclusion of the nozzle guide vane in both the HPG and LPA ducts, which help in maintaining the inlet velocity triangles and effectively preventing premature flow turning. These guide vanes have been neglected in the optimisation campaign in order to minimise computational effort.

Despite an increase in predicted shaft power, the optimised design still features large periods of positive power, such as between $0^\circ$ and $10^\circ$ azimuth angle. This implies further potential for optimisation. The presented optimisation study focused on the rotor shape using an existing port arrangement, port angles, and opening heights. It can therefore be deduced that the current optimisation is too constrained in terms of parameters to provide the best possible overall design.

The candidate design chosen from the pool of design iterations features a 3% lower pressure ratio in favour of greater power extraction. A decrease in pressure ratio hints towards an impeded ability to exchange energy between high- and low-pressure gas streams. To further investigate this, entropy generation rate has been integrated across the volume of a single channel using (16) and recorded as the channel rotates. The results of this are shown in Fig. 17a. The distribution shows that the largest amount of entropy is generated when the channels are exposed to the HPG duct. Subsequently, the entropy production reduces until the point, where the channels are exposed to the LPA duct, where low-pressure low-temperature gas is flowing into the channels. An additional contribution is given by channel-to-channel interaction. This accounts for a gradual increase in entropy production from $70^\circ$ azimuth angle onwards. Applying (13) and computing compression and expansion efficiency for both baseline and optimised model shows that the surplus in entropy is reflected in a reduced ability to compress and expand the gas streams efficiently through. This accounts for a 5% points decrease in the computed efficiencies, as shown in Fig. 17b. The same trends are witnessed when integrating entropy production rate over an entire cycle, albeit in a stronger fashion.

Fig. 18

Detailed view of instantaneous power output (top) and the associated incidence angle comparing baseline and power optimised design

The increase in entropy production is particularly prominent during the opening of the HPG duct. This is the phase with the largest instantaneous power generation, as shown in the top graph in Fig. 18. Evaluation of the corresponding incidence angle during this period is exhibited in bottom graphs in Fig. 18. The discrete manner in which the channels are exposed to the ducts results in a large variation of incidence angle. During the initial phase of channel exposure to the HPG duct at $\theta _1$, the data suggest that while the incidence angle distribution is similar between baseline and optimised design, the incidence angle for the baseline design reaches positive values earlier than the optimised design. Subsequently, the incidence angle evolves in the phase between $\theta _1$ to $\theta _2$ to moderately positive incidence angles of up to $3.8^\circ$. During this stage, a single channel commences to generate shaft power. It continues to do so even though the incidence angle then changes again, at $\theta _3$, to negative angles of $-\,8^\circ$ to $-\,10^\circ$. Despite the negative incidence angle, the channel continues to generate power, which is attributed to the non-axial camberline shape. Finally, a local peak in shaft power generation occurs as the HPG duct is in the process of closing and as incidence angle increases continuously.

Fig. 19

Contour plots comparing a static pressure distribution with the locations for primary and secondary shock waves; b mass fraction of gas entering through the HPG port and c normalised local entropy generation rate for the high-pressure zone of baseline (left-hand side) and optimised geometry (right-hand side)

Fig. 20

Comparison of shock patterns for baseline (left) and optimised (right) designs. Static pressure distribution along the channel length in the high-pressure zone at four different time instances (${\theta _1 {-} \theta _4}$) exposes implications of the changes in geometry on the generation of primary and secondary shock waves

Fig. 21

Comparison of a pressure distribution and b mass flow rates at the inlet port side (upper plot) and outlet port side (lower plot) for baseline (solid, black curve) and optimised geometry (broken, red curve)

The contour plots presented in Fig. 19 represent a single snapshot in time located in between channel hub and shroud. The given view focuses on the high pressure region in between the HPG and HPA duct comparing the baseline model (left-hand side) with the optimised geometry (right-hand side). Figure 19a displays static pressure for baseline and optimised geometry. Locations of the right-running primary shock waves are indicated (marker a1). In case of the optimised design shown on the right-hand side, it appears that the shock strength of the primary shock wave is slightly reduced compared to the baseline (left-hand-side plot). The secondary shock wave strength, on the other hand, remains unchanged (marker a2).

For the given loop flow ratio $\lambda =1.7$, the cold flow penetration length is below 50% of the channel width. This is marked (b1) in mass fraction contour plot in Fig. 19b. A mass fraction of 1 signifies flow entering the wave rotor through the HPG port, while a mass fraction of zero denotes flow entering through the LPA port. This facilitates the analysis of mixing and interaction of cold and hot streams within the channels. As $\lambda$ was kept constant, there is no change in the scavenging characteristics for the optimised geometry. In both geometries, the interface between hot and cold gases is highly skewed and three-dimensional (marker b2). In addition, the axial leakage path fosters mass transfer into the channels prior to the channels being exposed to the ports (marker b3).

Generation of irreversibilities shown in Fig. 19c is achieved by applying (14) and normalising the data for better visual representation for both baseline and optimised version by the same arbitrary number. Through comparison of the local entropy generation rate with the pressure and mass fraction field, it transpires that the bulk of the losses are generated at the demarcation surfaces between the hot and cold gas streams (marker c2). These losses appear to prevail over losses generated across the shock waves. Further losses are incurred as a consequence of turbulence and flow separation as the flow enters from the port into the channels (marker c1). Hence, the increased torque extraction of the optimised model appears to take place in conjunction with an increase in entropy production. Overall, the local results support the observations made in Fig. 17a, b attributing greater entropy production to the optimised rotor.

Figure 20 shows further insight into the change in gas dynamics between baseline and optimised geometry and the effect of the change in geometry and channel lengths on the shock patterns and pressure distributions. The left-hand-side contour plot represents the pressure distribution in the baseline geometry, while the right-hand-side contour plot denotes the optimised geometry. Qualitatively, the contour plots indicate differences, in particular visible in the high-pressure zone in the vicinity of the HPG and HPA ports.

Furthermore, extracting the static pressure distribution within the channel in the high-pressure zone (i.e., during opening of HPG and HPA) at four different instances $\theta _1{-}\theta _4$ reveals ramifications on the generation of primary and secondary shock waves. Data extracted from the baseline model (solid, black line) are plotted against data from the optimised geometry (broken, red line).

Fig. 22

Experimental results for a power and b pressure ratio comparing baseline and optimised rotor. The data are plotted against inlet mass flow rate through the HPG duct and were normalised by the baseline data at the target optimisation point

Fig. 23

Comparison of performance difference between the baseline model and optimised design for Q2D CFD and 3D CFD and as witnessed in the experiments. a depicts the power increase and b the absolute pressure ratio difference expressed in percentage

Initially, upon opening the channel to the HPG port (at $\theta _1$), increased incidence losses for the optimised geometry reduce the overall strength of the right-travelling primary shock wave manifesting itself in a 6% decrease in shock strength. Due to the 5% smaller channel length for the baseline channel shape, the primary shock wave reaches the outlet side sooner compared to the optimised design. The distribution in Fig. 20 at $\theta _2$ shows that for the baseline geometry the wave is reflected back as a secondary shock wave. In case of the optimised geometry and at the same instance in time, the wave is still travelling to the right just about to impinge on the outlet side. Finally, time instances $\theta _3$ and $\theta _4$ illustrate the generation of the reflected secondary shock wave. Again, as the secondary shock wave is generated earlier in case of the baseline model, there is a time shift between baseline and optimised model. Nonetheless, the overall shock strength in both examples is comparable and appears to be unimpeded by the difference in geometry.

The pressure distribution in the ports through a single cycle at both inlet and outlet sides is shown in Fig. 21. Despite the alterations in rotor geometry, there appears to be little qualitative difference between baseline and optimised geometry in terms of the pressure distributions in the inlet and outlet ports. Similarly, looking at instantaneous mass flow rates at inlet and outlet sides exhibits no significant changes in flow rates. The mass flow rates through the ports depict that the predominant flow direction is indeed in axial direction, with no significant instances of flow reversal. This situation persists for the optimised geometry.

Therefore, despite increased incidence and mixing losses determined from the loss audit, the modifications made to the channel shape appear not to have changed the overall gas dynamics significantly. The most prominent penalty can be witnessed with respect to the generation of the primary shock wave and altered impingement timing of the shock waves at the end plates. This is a consequence of larger incidence losses and is responsible for the decrease witnessed in pressure ratio. Nonetheless, secondary shock wave strength remains unaltered along with the pressure and mass flow distribution along the stator inlet and outlet sides.

Fig. 24

Inlet guide vanes located in the high- and low-pressure ports on the inlet side to maintain inlet velocity triangles

3.3 Experimental validation

Both rotor designs were tested experimentally on the University of Bath gas stand at the same operating conditions as imposed as boundary conditions in the CFD computations. Figure 22 shows a back-to-back comparison of the two wave rotor designs for power output and pressure ratio. At the target optimisation point of 32 g/s, the data indicate a substantial increase in shaft power output of 74% and a slight decrease in pressure ratio of approximately 4%. Figure 23 yields a comparison of the relative difference in performance between the baseline model and the optimised model for power (Fig. 23a) and pressure ratio (Fig. 23b) for the two numerical models used (Q2D CFD and 3D CFD) and data obtained from the laboratory tests. The data yield that the trends witnessed in the numerical models translate into experimental testing. The laboratory test data exhibit an increase in power by a factor of approximately 1.74 and hence lower than the figures predicted by CFD. This can most likely be attributed to manufacturing uncertainties with respect to the inlet guide vanes, as shown in Fig. 24. This exhibits deviations from the ideal and smooth shape used in the numerical campaign. In terms of pressure ratio, similar trends are visible. There clearly is a pressure ratio penalty for the optimised design, which increases from approximately 3% in the numerical model to 4.0% witnessed in the experiments. The reason for the increase is most likely found in the fact that the CFD model does not consider any interaction of the flow within the channel with the leakage cavity. Nonetheless, the experimental data give fidelity in both the Q2D and 3D model to be able to reproduce trends with reasonable accuracy without the need to resort to higher-fidelity models and greater domain detail.

4 Conclusions and outlook

In this paper, a numerical optimisation of an axial throughflow wave rotor turbine is presented for the first time. The approach employed a combined quasi-two-dimensional, transient CFD, and hybrid evolutionary algorithm. The primary objective of the optimisation campaign was to improve the power output while ensuring pressure exchange capabilities are not significantly penalised.

A hybrid GA has been proposed that first runs a predefined number of generations using a standard elitist GA. The generated data are then used to train a surrogate model, which is then used to perform auxiliary optimisation. The surrogate model is then constantly updated as the optimisation routine minimises it using additional constraints. These constraints select the subsequent design candidate to undergo costly CFD simulations based on its distance from previously evaluated points. The optimisation then runs through a defined search pattern: varying exploitation (local search) and exploration (global search).

For the wave rotor optimisation, a single operating point was selected and the rotor channel camberline shape and wall thickness parameterised using five different parameters. The stator port solution and port angles were kept constant. In addition, the number of channels in the rotor was allowed to vary as well. The objective function aims at increasing power output with additional constraints on minimum pressure ratio, channel width, and the number of channels.

Two optimisation runs were started, where one employed the hybrid-GA algorithm and the other MATSuMoTo. Both methods converge quickly to a minimum with the hybrid GA requiring fewer than 100 design evaluations and MATSuMoTo approximately 350. The best candidate design of both optimisations indicated an increase in predicted power output by a factor of 1.78.

In order to evaluate the designs further, a direct comparison of the optimised design with the baseline design in three-dimensional CFD confirmed a higher channel loading and thus an increase in power output. The three-dimensional CFD results confirm the findings of the quasi-two-dimensional model emphasising the validity of using such a reduced-order model for optimisation.

Table 8 List of mathematical functions used as test problems to compare a standard GA, PSO, MATSuMoTo, and the hybrid algorithm

The results imply that the increase in momentum transfer results in an increase in entropy production, in particular at the HPG port inlet and at the highly skewed mixing interface between hot and cold gas streams. In addition, the evolution of inlet incidence angle exhibits that the incidence angle for the optimised design is initially more negative during the opening phase of the channel to the HPG duct. This accounts for a higher flow separation and thus entropy production. Overall, incidence angles appear to vary substantially as a consequence of finite passage opening and are predominantly negative throughout the phases with maximum shaft power generation implying that camberline curvature is the dominant driver to shaft power generation as opposed to incidence momentum change. However, the results indicate that further optimisation potential may be exploitable if the stator arrangement is included in the optimisation.

A gas dynamic analysis within the channels revealed that the geometry modification affects the generation of the primary shock wave as well as the timing at which waves impinge at inlet and outlet sides. The effect on primary shock wave generation is the main driver for the witnessed decrease in pressure ratio, while secondary shock waves and port pressure and mass flow distributions remain unchanged by the optimisation results.

Comparing baseline and optimised design on a gas stand indicates an increase in power output of approximately 1.74 and decrease in pressure ratio of 4% confirming the trends witnessed in the numerical models. It can thus be concluded that the quasi-two-dimensional represents a viable tool to efficiently optimise wave rotor designs.

Appendices

Appendix 1

1.1 Numerical test cases

Table 9 Average number of function evaluations necessary to be within 1% of the global optimum for different optimisation algorithms on the test problems

In order to obtain an estimate of the performance and gain fidelity in the routine, the hybrid method is tested on a number of synthetic optimisation problems for global optimisation. This includes the Dixon and Szegö test functions [60] (Branin, Goldstein-Price, Hartman 3, Hartman 6, Shekel 5, Shekel 7, and Shekel 10) as well as the non-convex nonlinear Rosenbrock, Ackley, and Rastrigin functions. A summary of the test problems is given in Table 8. Furthermore, the performance of the hybrid algorithm is compared with a standard elitist GA, PSO, as well as MATSuMoTo, another surrogate-based optimisation routine [61,62,63,64]. MATSuMoTo differs as it uses design of experiments (such as Latin-hypercube sampling) to create the initial training set and uses radial basis functions (RBF) to construct the surrogate model. Instead of a user-defined search pattern, MATSuMoTo uniformly selects data points from the design space and perturbing them with a randomly chosen standard deviation. For each of the candidate points, a weighted score based on their objective function value predicted by the surrogate model is computed and the candidate with the best score is selected as the subsequent evaluation point. Furthermore, the algorithm cycles through a predefined weight pattern in order to alternate between local exploitation and exploration.

The idea behind using the test problems is that they represent complex surfaces that pose a reasonable challenge for any optimisation routine with non-convex shapes, multiple parameters, and multiple local minima. They may not be able to recreate the challenges for wave rotor optimisation, but are cheap to solve and should at least give an idea whether an optimisation routine has the potential for swiftly finding optimal solutions. Therefore, each optimisation routine is run through the test problems 30 times with a fixed computational budget of 1000 expensive function evaluations. The global optimum as well as its location of each function is known, so that the optimisation can be regarded as successful when the best confirmed function value falls within a threshold of 1% of the global optimum. In case the function value at the global optimum is zero, the threshold is set to 0.05. For the hybrid algorithm, three different search patterns are examined with the aim of finding a universally applicable pattern that can be further used for optimising the wave rotor turbine. The search patterns used feature two of length six and one of length five $\beta _1=[0.9, 0.75, 0.25, 0.05, 0.03, 0]$, $\beta _2=[0.95, 0.25, 0.05, 0.03, 0]$, and $\beta _3=[0.95, 0.5, 0.25, 0.1, 0.05, 0]$.

Fig. 25

Box plots showing the statistical data for the number of function evaluations to reach the global optimum with a relative error of 1%. For runs with failed cases, no box plot is provided. a refers to the results for the two-dimensional Branin function, b to the two-dimensional Goldstein-Price function, c the three-dimensional Hartman function, d to the two-dimensional Rastrigin function, e to the six-dimensional Hartman function, and f to the four-dimensional Rastrigin function

Fig. 26

Comparison of the convergence rate of the hybrid GA and MATSuMoTo with other hybrid optimisation routines for a the five-dimensional Rastrigin function and b the two-dimensional Rosenbrock function [65, 66]

For each test problem, GA parameters of population size and crossover fraction were varied in order to estimate the sensitivity of results to the choice of these parameters. The same approach was done for PSO regarding swarm size, self-adjustment, and social adjustment weights. In both instances, only the results of the best set of options are used for comparison of the different optimisation methods. The results for these cases are reported in Sect. 4.

Appendix 2

1.1 Results of numerical test cases

The average number of function evaluations to reach a relative error of 1% of the known global minimum for GA, PSO, hybrid algorithm, and MATSuMoTo is given in Table 9. In addition, the standard deviation is presented in parentheses along with the number of instances, when the algorithm was unable to detect the global optimum within 1000 iterations. In case a certain number of runs do not finish within the iteration limit, the run is excluded from being evaluated in the averaging process. For all pure GA runs, a population size of ten was used and a crossover fraction of 0.6 was used, with the exception of Shekel 5, which required a crossover fraction of 0.5. Pertaining to PSO, the general trend across all trials implied lower swarm sizes leading to accelerated convergence. Therefore, the swarm size chosen varied in between 5 and 10 depending on the problem in conjunction with self-adjustment and social adjustment weights of varying between 1.49 and 1.59, respectively.

Nonetheless, even with a tweaked set of model parameters, the elitist GA and PSO, as implemented in MATLAB R2017b, struggle to consistently detect the optimum over 1000 function evaluations. The two evolutionary algorithms face particular challenges with multidimensional functions containing a large number of local minima. For such problems, such as higher-order Rastrigin functions, higher numbers of failed searches are recorded. Overall, PSO appears to perform better with fewer failed cases and faster detection.

All hybrid variants as well as MATSuMoTo outperform GA and PSO. Among the different search patterns used for the hybrid GA, $\beta _2$ and $\beta _3$ perform better on all test functions with the exception of the two parameter Rastrigin function. Search pattern $\beta _2$ performs best on the two-dimensional Branin, Goldstein-Price, and Ackley functions, but exhibits comparably poor results for the multidimensional Rastrigin function as well as the Rosenbrock test. In general, it appears that search pattern $\beta _3$ performs marginally better than either of the two other search strategies, in particular for higher-order multimodal functions such as the four-dimensional Rastrigin function.

Compared with the surrogate-based optimisation of Müller’s MATSuMoTo using thin plate splines RBF, the hybrid approach performs better on the Branin, Goldstein-Price, Shekel 5, three- and four-dimensional Rastrigin as well as the Ackley function. It is particularly surprising that MATSuMoTo, at least with the proposed set-up, struggles with the higher-order Rastrigin problems resulting in eleven failed cases. Both approaches, however, yield similar results for Hartman 3, Shekel 7, and Shekel 10. Applied to the multidimensional Hartman 6, two-dimensional Rastrigin, and the Rosenbrock functions, MATSuMoTo performs better.

Figure 25 shows box plots for a selection of test problems. Only cases without any failed trails are shown. Examining the data for Branin, Goldstein-Price, Hartman 3, and the two-dimensional Rastrigin demonstrates a wide spread of between the $25{\mathrm {th}}$ and $75{\mathrm {th}}$ percentile for both GA and PSO with a number of outliers in particular seen for the GA applied to the Branin problem. Overall, the data show that PSO may converge faster than the standard elitist GA. Comparing the hybrid method with MATSuMoTo, on the other hand, proves to achieve convergence with little variation. However, there is still the possibility of outliers requiring a higher number of iterations for both surrogate-based models.

Finally, to see how the hybrid GA and MATSuMoTo fare against other recently introduced hybrid methods, two single trials are run on the Rastrigin function with five parameters and the two-dimensional Rosenbrock function. The comparison is made with another evolutionary algorithm, namely a classical differential evolution (DE) algorithm, hybridised with a steepest descent method developed at the von Karman Institute [65, 66]. The convergence rate for both test problems is presented in Fig. 26 showing the best detected function value for Rastrigin in (a) and Rosenbrock in (b). Although the hybrid DE versions converge steadily and quicker compared to the original DE, both search patterns of the hybrid GA exhibit steeper convergence and being able to detect the global optimum after approximately 1200 iterations. For the Rosenbrock function, hybridised DE again shows a steady convergence rate, while the hybrid GA with search pattern $\beta _3$ quickly determines the global optimum after less than 60 iterations. Both hybrid GA with search pattern $\beta _2$ and MATSuMoTo exhibit extended constant periods, with the former finally reaching the global minimum after around 1000 iterations. The latter, however, does not manage to locate the optimum within the imposed 1500 iterations in the particular shown run.

In sum, the statistical evaluation implies that combining evolutionary algorithms with surrogate models or substituting evolutionary algorithms completely can yield significant performance gains making these approaches particularly attractive for expensive CFD simulations, such as required for wave rotors.