**CFD Optimization Process of a Lateral Inlet**/**Outlet Di**ff**usion Part of a Pumped Hydroelectric Storage Based on Optimal Surrogate Models**

#### **Xueping Gao 1, Hongtao Zhu 1, Han Zhang 1, Bowen Sun 1,\*, Zixue Qin <sup>1</sup> and Ye Tian <sup>2</sup>**


Received: 5 March 2019; Accepted: 4 April 2019; Published: 10 April 2019

**Abstract:** The lateral inlet/outlet plays a critical role in the connecting tunnels of a water delivery system in a pumped hydroelectric storage (PHES). Therefore, the shape of the inlet/outlet was improved through computational fluid dynamics (CFD) optimization based on optimal surrogate models. The CFD method applied in this paper was validated by a physical experiment that was carefully designed to meet bidirectional flow requirements. To determine a good compromise between the generation and pump mode, reasonable weights were defined to better evaluate the overall performance. In order to find suitable surrogate models to improve the optimization process, the best suited surrogate models were identified by an optimal model selection method. The optimal configurations of the surrogate model for the head loss and the velocity distribution coefficient were the Kriging model with a Gaussian kernel and the Kriging model with an Exponential kernel, respectively. Finally, a multi-objective surrogate-based optimization method was used to determine the optimum design. The overall head loss coefficient and velocity distribution coefficients were 0.248 and 1.416. Compared with the original shape, the coefficients decrease by 6.42% and 40.28%, respectively. The methods and findings of this work may provide practical guidelines for designers and researchers.

**Keywords:** pumped hydroelectric storage; inlet/outlet; surrogate model selection; multi-objective optimization process

#### **1. Introduction**

In recent decades, many electricity generation problems have been caused by their high production cost and environmental impact [1]. It is generally accepted that renewable energy can provide promising solutions for the problems, and much research has been focused on numerous techniques, such as hydropower [2–4], fuel cells [5], biofuels [6], wind power [7–9], and solar power generation [10]. Pumped hydroelectric storage (PHES) has attracted widespread interest because of its great flexibility and storage capacity in improving grid stability of other renewable energy sources [11,12]. The lateral inlet/outlet plays a critical role in the connecting tunnels of the water delivery system in a PHES, and its flow conditions can affect the economic benefit and operation stability of the PHES to a great extent. The PHES operates under pumped or generation conditions according to the command, and the water flows between the upper and lower reservoir through the inlet/outlet. Compared with ordinary hydropower stations, the hydraulic requirements for the inlet/outlet structure of the PHES require more attention.

Much research in recent years has focused on the hydraulic performances of the inlet/outlet experimentally and numerically. Müller et al. [13] carried out prototype measurements by Acoustic Doppler Current Profiler (ADCP) to study flow velocities in the reservoir; subsequently a numerical model was built to further understand the flow development near the inlet/outlet. In order to ensure satisfactory hydraulic performances, Bermúdez et al. [14] conducted numerical and physical model studies on the inlet/outlet of the Belesar-III power station in Northwest Spain. The initial design was optimized by numerical method, and the new design presented a more homogenous flow distribution and a reduced head loss. Cai et al. [15] performed an experimental study on the effects of the shape and position of separation piers on the flow characteristics with two operations. Gao et al. [16] tried to improve the velocity distribution at the intake-outlet orifices through 20 types of shape optimization experiments. Sun et al. [17] found the incorrect arrangement of separation piers probably contributes to disadvantageous velocity distribution for the trash rack. Ye et al. [18] simulated turbulent flow in a lateral inlet/outlet to explore the diffusion segment shape that can obtain better velocity and discharge distributions. Based on the validation of the experiment, flows in the lateral inlet/outlet were numerically simulated focusing on the variations in discharge distributions over different orifices; the results showed that the flow non-uniformity was improved successfully by adjusting the diffusion segment and width between the separation piers [19]. The investigations indicate that computational fluid dynamics (CFD) has been used frequently to optimize the structure parameters of the inlet/outlet that can easily influence the hydraulic objectives. However, the process mentioned above is a trial-and-error approach requiring intensive work and a lot of time [20], and the results may be subjective to some extent.

CFD coupled with optimization techniques has been used to determine the best design scheme for wind turbine devices [21–23], and the adopted methods provide a reference for inlet/outlet optimization. For complex optimization problems, two methods have been generally applied to achieve optimization. When there are many individual objective functions, the functions can be combined into a single composite one; nonetheless it can be difficult to determine the proper and typical selection of the weights. In addition, the combination method would provide only a single result that cannot be chosen to make a trade-off. Therefore, multi-objective evolutionary algorithms (MOEAs) are more desired by decision makers. For MOEAs, it is common to use the concept of Pareto solutions. Pareto solutions consist of many non-dominated solutions where some gain in objectives always causes some sacrifice to others. MOEAs, such as Non-dominated Sorting Genetic Algorithm II (NSGA-II) [24,25], can offer many Pareto results in one single simulation process and have been extensively applied in engineering optimization.

Although MOEAs can obtain many solutions, these algorithms also require an even larger number of simulations and an unacceptable computation time in practice. Thus, multi-objective optimization based on the surrogate methods is receiving increasingly more support and attention. These approaches present objective functions by surrogate models, replacing expensive high-fidelity design simulations. The model is an effective engineering technique that can considerably save computing expense while still providing reasonably accuracy. Therefore, the multi-objective surrogate-based optimization method makes it more feasible to perform exhaustive global searches in the design space for complicated shape optimization. The optimization result will be more satisfactory because the exploration is not subjected to a limited number of simulating data and computing resources. General surrogate models contain Response Surface Methodology (RSM) [26,27], Radial Basis Function (RBF) [28–32], and Kriging [33–35]. However, owing to the diverse functional characteristics of the popular models, it is important to choose the most appropriate models for the responses of hydraulic objectives to diffusion segment parameters. For that reason, a recently surrogated model selection method proposed by Mehmani et al. [36], called the Concurrent Surrogate Model Selection (COSMOS), was adopted to construct suitable models in this paper.

Although a multi-objective CFD optimization process for the inlet/outlet was developed in previous works [37,38], there are some drawbacks in the original process: (1) The CFD method was indirectly validated only by comparing the result with that of a published reference [18], which weakened the reliability of the optimization process. (2) Only a single surrogated model (RSM) was considered to accelerate the optimization procedure. However, there is no universal surrogate model for every practical problem, and suitable model selection is critical for accuracy. (3) The objectives for the pump and generation mode were simply added up to estimate the total performance. Nevertheless, the significance of the pump mode and generation mode is probably not equal in many cases. Therefore, in this study a physical model experiment was well-designed to meet the requirement of bi-directional flow, and a more extensive validation of the CFD model was performed. Based on the validated CFD model, the COSMOS method was applied to select suitable surrogate models for the objectives of the inlet/outlet, to further improve the accuracy of the original optimization process. In addition, in order to better evaluate the overall hydraulic performance, the selection of objective functions was reconsidered and weight coefficients of different operation modes were introduced. In the present article, the shape of an inlet/outlet was improved through CFD optimization based on optimal surrogate models. The aim of this work was to apply the optimal surrogate optimization process to determine the optimum shape of the inlet/outlet, to simultaneously obtain better hydraulic performances under bidirectional flow conditions.

#### **2. Problem Description**

The typical shape of the lateral intake is shown in Figure 1. When the PHES is in the generation mode, water flows into the lower reservoir through the tunnel and the inlet/outlet structure. Because water flows out from the inlet/outlet in this mode, the phenomenon is also called "outflow". Along the outflow direction, the structure generally contains a transition part (a structure connecting a circular cross-section tunnel and a diffusion part), a diffusion part, a rectification part, and an anti-swirl part.

**Figure 1.** General structure of the lateral intake-outlet.

The diffusion part is very sensitive to the performance of the intake, and it is generally improved for the structure design. Figure 2 present flow passages and basic parameters in the diffusion part. This parameters include the length of the diffusion part (*L*DS),height of the inlet/outlet orifices (*H*I/O),wide of the separation pier in the middle orifices (*W*DSM),width of the separation pier in the side orifices (*W*DSS),horizontal diffusion angle (η), and vertical angle (θ). The left side displays the 3D model for the flow passage. On the right side, plane sections of the model are demonstrated, and the flow passage and separation piers are marked out with light grey and dark grey, respectively. It should be noted that the tunnel diameter D in the PHES is generally decided beforehand by the difference of the water head between the upper and lower reservoir, the installed capacity of turbines, topographic and geological conditions etc. For the inlet/outlet optimization, the value of D is fixed (D = 7.2 m in the current study) because a change may need adjustment of the whole layout of the water conveyance system resulting in excessive costs. Besides, the number of orifices (4) is also fixed in view of it being the most commonly used form of the lateral inlet/outlet in the PHES of China. For more information about the shape parameters, readers can refer to the previous work [38].

**Figure 2.** Basic parameters for the diffusion part.

The lateral inlet/outlet hydraulic performance is characterized by three parameters: the head loss, the velocity, and the flowrate distribution [37–39]. According to the Bernoulli equation, the head loss is computed as:

$$h\_{0\cdot 1} = \nabla\_0 - \nabla\_1 - \frac{av^2}{2g} \tag{1}$$

$$h\_{1\cdot 0} = \nabla\_1 - \nabla\_0 + \frac{av^2}{2g} \tag{2}$$

$$
\xi = \frac{2gh\_f}{av^2} \tag{3}
$$

where the head loss *hf* = *h*0-1 and *hf* = *h*1-0 are used for the inflow and outflow mode, respectively; ξ is the coefficient of the head loss; ∇<sup>0</sup> and ∇<sup>1</sup> are the piezometric head in the reservoir and tunnel; *v* is the bulk velocity at the boundary of the tunnel; α is the kinetic energy correction coefficient.

Numerous vibration-caused trash rack failures have occurred because the trash rack may be exposed to high velocities at the orifice cross-section [40]. Hence *C*<sup>V</sup> is proposed as a typical index evaluating the uniformity of the flow velocity at orifices [39,41], and a lower *C*<sup>V</sup> implies a better velocity distribution in the inlet/outlet. Considering there are four orifices divided by the separation pier, *C*<sup>V</sup> could change with different orifices due to the influence of horizontal diffusion. In order to represent the overall performance of the velocity distribution, the worst case i.e., the maximum value of *C*<sup>V</sup> among the four orifices, is selected. Thus *C*<sup>V</sup> is defined as:

$$C\_V = \max(\frac{v\_{\max,i}}{v\_i})\tag{4}$$

where *v*max is the maximum value among the velocities at the orifice cross-section; *i* indicates the index of orifices; *v* indicates the average velocity at the orifice cross-section.

In order to further improve the discharge distribution in the horizontal direction, the separation piers are designed to divide the flow passage into sub-tunnels. Then *C*<sup>Q</sup> is used to estimate the degree of flux uniformity by comparing the discharge between adjacent orifices, which can be defined as:

$$C\_Q = \frac{\max(Q\_{i,i+1}) - \min(Q\_{i,i+1})}{\min(Q\_{i,i+1})} \times 100\% \tag{5}$$

where *Qi* indicates the rate of flow at the different orifices (*i* = 1, 2, 3).

For the purpose of determining a good compromise between the two modes of operation, these parameters under inflow and outflow conditions should be carefully combined. In a recent study, Gao et al. [38] simply added up the parameters for the pump mode and the generation mode to estimate the total performance. However, the performance of the generation mode is probably more significant than that of the pump mode in many cases. The pumped water can be stored in the upper reservoir under the pump mode during the low load period of the power grid, while under the generation mode the water will flow into the lower reservoir to generate electricity during the peak load period. Designers tend to pay more attention to the head loss under the generation mode in order to achieve higher economic benefits. Therefore, a more reasonable approach is to use weights to deal with the combinations of the two modes. For that reason, the new combined objectives to evaluate the overall performance of the inlet/outlet are defined as follows:

$$
\omega\_{\text{total}}^{\mathbb{K}} = \omega\_1 \xi\_{\text{in}} + \omega\_2 \xi\_{\text{out}} \tag{6}
$$

$$\mathcal{L}\_{\rm V,total} = a\nu\_1 \mathcal{C}\_{\rm V,in} + a\nu\_2 \mathcal{C}\_{\rm V,out} \tag{7}$$

where ξtotal and *C*V,total are the overall head loss coefficient and velocity distribution coefficient, respectively; ω<sup>1</sup> and ω<sup>2</sup> are the weight coefficients of the objective function under the conditions of inflow and outflow, and ω<sup>1</sup> + ω2= 1. In this article, the values of ω<sup>1</sup> and ω<sup>2</sup> are taken as 0.33 and 0.67 based on the specific engineering designer's suggestion. However, they can be flexibly selected according to different projects.

Another drawback of their study [37,38] is the definition for the overall *C*Q, which is expressed as the difference of *C*<sup>Q</sup> between the pump and generation mode. *C*Q,in is the discharge uneven distribution coefficient under inflow condition, *C*Q,out is the discharge uneven distribution coefficient under outflow condition. When it is regarded as an objective function, it seeks to minimize the difference between *C*Q,in and *C*Q,out, rather than *C*Q,in and *C*Q,out themselves. If there are some bad points in the design space, on which the values of both *C*Q,in and *C*Q,out may be large, the difference (i.e., the overall *C*Q) may be small. Besides, according to the guideline the discharge distribution is sufficiently uniform when *C*<sup>Q</sup> is less than 10% [39]. As a result it is more appropriate to choose *C*<sup>Q</sup> as the constraint, and the head loss and velocity distribution (ξtotal, *C*V,total) can be optimized more flexibly.

Therefore, the objective functions were set to find the minimum values for ξ and *C*V, while *C*Q, η, and θ were treated as constraint conditions. Table 1 lists the design parameters and constraints in this study.

**Table 1.** Design parameters and constraints.


#### **3. Numerical Model and Validation**

#### *3.1. Numerical Model*

Figure 3 shows a detailed view of the structured grids in the research area. To increase the accuracy of the simulation, boundary layer grids were added to ensure that the dimensionless wall distance *y*+ of the first points near the walls was between 30 and 300. The distribution of *y*+ values for the inlet/outlet is shown in Figure 4.

**Figure 3.** Grids generation for the computational domain: (**a**) grids on the boundary; (**b**) grids at the orifice cross-section, and (**c**) grids on the orifice horizontal section.

**Figure 4.** The Distribution of y+ values for the inlet/outlet: (**a**) the inflow condition; (**b**) the outflow condition.

Considering the accuracy and efficiency of the CFD simulations, it is important to work out a reasonable number of mesh elements [42–44]. Based on a previous grid independence study [38], a grid amount of 1.4 <sup>×</sup> <sup>10</sup><sup>6</sup> was sufficient for CFD simulation, and similar grid sizes were thus adopted in this article.

As presented in Figure 1, a three-dimensional computational domain far larger than the inlet/outlet structure was developed to accurately simulate the flow pattern in the inlet/outlet. As to the boundary conditions, a uniform velocity was prescribed at the tunnel boundary, and a pressure condition was imposed at the reservoir boundary. The flow discharge changed slightly according to the operating conditions of the reversible turbines. A summary of boundary conditions for the Computational Fluid Dynamics (CFD) case under the generation and pump mode is listed in Table 2.


**Table 2.** Boundary conditions applied in the inlet/outlet CFD simulation.

In the present study, the Reynolds Averaged Navier Stokes (RANS) equations were solved using ANSYS Fluent 17.0 [45]. The simulation was performed in a segregated manner adopting the realizable *k-*ε model combined with the wall function. The Navier–Stokes equations were solved by means of pressure implicit with splitting of the operator's algorithm, and a second-order upwind method was employed to realize the discretization. The converged solution was obtained after the residual levels became less than 1 <sup>×</sup> 10−4. For detailed information about the governing equation and numerical settings, the reader can refer to references [38,45,46].

#### *3.2. Experimental Setup*

The CFD method in the original study [38] was validated by comparing the result with that of a published reference [18]. Although the results were found to be in good agreement, there were obvious differences in the velocity distribution between the case in the reference [18] and the case in this study, which weakened the reliability of the validation in the original study. Therefore, the CFD method still needs to be validated further, and a physical model experiment adopting a scale of 1/60 was conducted in this paper. As shown in Figure 5, an experimental setup was established which was 5.2 m long, 1.2 m wide, and 0.6 m high.

**Figure 5.** The experimental setup of the inlet/outlet.

In order to show the experiment more clearly, the side view of the layout plan is illustrated in Figure 6. Control valves and electromagnetic flowmeters with an accuracy of 0.5% for the measured ranges were used to control the flow discharge. By switching the different control valves, bidirectional flow conditions could be realized in the experiment. In the inflow condition, the inflow control valves were open, while the outflow control valves were closed. Water from the high constant-head tower flowed into the reservoir channel through the pipe. Before the water entered the channel, a device was set up to align the flow pattern. Finally water flowed into a lower collecting tank, from which it was pumped back to the constant-head tank using a centrifugal pump. While in the outflow condition, the inflow control valves were closed, and the outflow control valves were open. Water was conveyed to

the reservoir tank through a horizontal pipe with a length 50D that was long enough to guarantee full development of the flow in the pipe. The discharge in the inflow condition was 4.359 L/s, corresponding to a Reynolds number Re = 46,277 in the tunnel; while in the outflow condition, the discharge was 5.795 L/s and the Reynolds number was 61,519.

**Figure 6.** The side layout plan for the experimental system of the inlet/outlet.

A Vectrino ADV instrument, developed by Nortek AS, was used to instantaneously measure stream-wise velocity of the diversion orifices at various depths. Based on the Doppler shift principle, the ADV instrument can measure velocities with ±1% accuracy in a measurement range of 1 mm/s. Each measurement point with a sampling volume of 7 mm worked at a sampling frequency of 50 Hz for 60 s of sampling time.

#### *3.3. Validation of the CFD Model*

The head loss coefficient in the inflow condition is obviously smaller than that in the outflow condition. In the inflow condition, the coefficients obtained from numerical simulation and experiment were 0.164 and 0.171 respectively, and their relative error was 4.09%. The coefficients of the simulation and experiment in the outflow condition were 0.315 and 0.338, and the relative error was 6.81%. Comparing the simulation results with the experimental data, although the numerical model slightly underestimates the coefficients, it is still acceptable to predict the results of the head losses.

For the purpose of comprehensively comparing the velocity distribution at the orifice cross-section (x = 10 m), there are three measured lines with equal distance at each orifice. As shown in Figure 7, the velocity data are obtained along the measured lines (A—A, B—B and C—C) at the cross-section.


**Figure 7.** Measured lines at the cross-section (x = 10 m).

Figures 8 and 9 compare the numerical solutions with the experimental ones for the velocity distribution in both conditions. In the legend, M and S indicate the middle and side orifice, respectively; A, B, and C indicate the measured lines A—A, B—B, and C—C, respectively. For example, M-A represents the measured line A—A at the middle orifice. The *y*-axis indicates the relative height of the

measured point normalized by the orifice height H, and the *x*-axis indicates the relative velocity also normalized by the average velocity of the measured line.

Figure 8a demonstrates that a good matching between the two approaches is obtained in the middle orifice. The results of both methods show that velocities among the three lines are very close to each other, and the distribution is also uniform along the height. In the side orifice, the numerical results are also in good agreement with the experimental ones. There are only slight differences in the velocity distribution among the measured lines. Compared with the flow pattern in the middle orifice, the more obvious change is that the mainstream appears at the lower part of the orifice, and the range of *y*/H is about 0.13–0.52.

While in the outflow condition, differences of the velocity distribution among the measured lines are also small in the side orifice. However, in the middle orifice the experimental results show that there are obvious differences among the velocities of three measured lines, which is also predicted by the simulation. In addition, the velocity distribution along the height is also not uniform, and negative velocities are generated at the top. For the M-B line, the *y*/H range of the negative region obtained from the experiment is about 0.8–1.0, and the numerical range is about 0.7–1.0. This difference is partly due to errors caused by the limited number of measured points in the experiment, and partly due to the slight overestimation for the strength of flow separation by the simulation. At the bottom of the orifice, the mainstream in the experiment is located a little higher than that in the simulation. Although there are some minor deviations in the outflow velocity distribution, the agreement is quite acceptable considering that the more complex flow separation occurs in the outflow condition.

**Figure 8.** Validation for the velocity distribution in the inflow condition: (**a**) middle orifice; (**b**) side orifice.

**Figure 9.** Validation for the velocity distribution in the outflow condition: (**a**) middle orifice; (**b**) side orifice.

#### **4. Optimization Methodology**

#### *4.1. Surrogate Models Selection*

The optimization algorithm coupled with the computationally intensive simulation requires intensive work and a lot of time. It is necessary to adopt the surrogate model because the CFD-optimization procedure requires too many computing resources, which weakens its application in practice. Therefore, surrogate models are employed to engage the optimization procedure in the paper. There are several surrogate model types as listed in Table 3, such as Response Surface Methodology (RSM), Radial Basis Function (RBF), and Kriging, which have been extensively used in recent years to solve engineering optimization design. Unlike other approaches, RSM applies simple algebraic expressions to generate the function relationship, and it is very suitable for complex engineering applications [47]. In addition, RBF has the best approach for high dimensionality optimization problems according to a literature survey [48]. For low dimensional optimization problems, the Kriging method can achieve good results in both accuracy and time. According to a survey of popular surrogate models in similar situations, RSM, RBF, and Kriging are considered as surrogate candidates in the researches.


**Table 3.** A list of shape optimization investigated using surrogate models in recent years.

When employing the surrogate model in the hydraulic optimization procedure, the challenge is how to create a model which is as accurate as possible. As mentioned above, RSM, RBF, and Kriging can be used to generate surrogate models, and their parameters are listed in Table 4. The kernel function and hyper-parameter also need to be chosen carefully, which requires an effective surrogate model selection method.

As we know, existing surrogate model selection methods generally contain three levels: selecting a model type, selecting a kernel function, and optimizing the hyper-parameters. Some researchers have used different error measures to separately select the model type, kernel function, or hyper-parameter [56–58]. With the intention of performing thorough selection, the Predictive Estimation of Model Fidelity (PEMF) was suggested by Mehmani et al. [59]. Furthermore, Mehmani et al. [36] advanced the Concurrent Surrogate Model Selection (COSMOS) on the foundation of the PEMF, and applied it to select suitable models by minimizing the model error. The error selection criterion is based on different situations and preferences. For the exploration of design space and parameter analysis, the median error can be used. On the other hand, it is more appropriate to choose the conservative maximum error for the structural safety optimization (e.g., to model vehicle crash simulation). Therefore, these criteria might be mutually conflicting or mutually promoting [36]. Considering the optimization for the hydraulic performance of the inlet/outlet may have certain requirements for the above two situations, the median and maximum errors of the surrogate model are treated as two selection criteria estimated by the PEMF method. Moreover, the framework will solve the multi-objective problem (minimizing the errors) to find the optimal models through the One-Step technique. In the technique, a single mixed integer nonlinear programming problem

(MINLP) can be solved by the NSGA-II. Considering different models have different numbers of kernels and hyper-parameters, the candidates are classed as three smaller categories according to the hyper-parameters [36]. The computation of the Kriging model is performed with the MATLAB toolbox DACE [60], as it has been widely used in many applications concerning the Kriging model. The problem is usually anisotropic, which is accounted for in the construction of the surrogate models.


**Table 4.** Parameters of surrogate model candidates.

#### *4.2. Optimization Algorithm*

Extensive literature based on the NSGA-II [61] has been reported in engineering applications (listed in Table 3), which is also chosen for this investigation. In NSGA-II, a random parent population is initially produced, and individuals are sorted according to rank and crowding distance. Then the optimum resolutions are chosen to generate offspring populations using genetic operators that are the same in the standard GA. Specifically, a binary tournament selection, Simulated Binary Crossover (SBX) and polynomial mutation are employed in NSGA-II. SBX simulates the binary crossover observed in nature, and the crossover distribution index determines how far away the children solutions are from their parents [62]. Deb et al. [63] suggested a polynomial mutation operator with a mutation distribution index that can control the amount of perturbation in a variable. Table 5 lists the parameters for the optimization process.

**Table 5.** The parameters for the optimization process.


To pick the optimum point from the Pareto frontier, Technique for Order Preference by Similarity to Ideal Solution (TOPSIS) is implemented in this article. TOPSIS, proposed by Hwang et al. [64], is a ranking method that can rank possible solutions on the Pareto frontier and select a trade-off optimum by measuring Euclidean distances.

#### *4.3. Optimization Process Scheme*

The multi-objective optimization process is built by combining the CFD numerical simulation, the design of experiment (DOE), the surrogate model selection based on the COSMOS, and the optimization method NSGA-II. The flow diagram of the process is presented in Figure 10.

First, the Optimized Latin Hypercube Sampling (OLHS) algorithm is applied in the Design of Experiment (DOE) to choose representative points. Then sample points are simulated using the CFD method, and the results are exacted to evaluate the objectives of the hydraulic performance. Based on the database of the DOE, a suitable surrogate model is constructed using COSMOS to improve the optimization efficiency and precision. The different surrogate models, such as RSM, RBF, and Kriging, are generated until the PEMF is deemed acceptable. Finally, NSGA-II is adopted to conduct the optimization. If the objective function does not change, or the maximum iteration number is reached, the optimization is terminated. Otherwise, the process will continue to loop with a new set of proposed design parameters.

**Figure 10.** Flow diagram of the multi-objective optimization process.

#### **5. Results and Discussion**

#### *5.1. Optimal Surrogate Models*

In the present study, COSMOS is applied to choose the appropriate surrogate models. First of all, OLHS as an efficient sampling strategy is applied in DOE sampling, in which 60 points are randomly selected to construct these surrogate models. The Predictive Estimation of Model Fidelity (PEMF) method is chosen to check the accuracies of these surrogate models.

By solving the MINLP problem with NSGA-II, the final results of the models for ξ, *C*V, *C*Q,in and *C*Q,out are presented in Figure 11. Φ*<sup>i</sup>* indicates the surrogate models containing *i* hyper-parameter(s). βH, βL, and β<sup>W</sup> indicate the correlation parameter of the design variable *H*I/O, *L*DS, and *W*DSM, respectively. A Pareto frontier is provided in each figure, and a trade-off solution is selected under comprehensive consideration of median and maximum errors. In Figure 11a, optimal parameters of the surrogate model of ξ are: the Kriging with a Gaussian kernel and β<sup>H</sup> = 2.20, β<sup>L</sup> = 9.41, β<sup>W</sup> = 0.11. As shown in Figure 11b, the optimal configuration of *C*<sup>V</sup> is: the Kriging model with an Exponential kernel and β<sup>H</sup> = 1.05, β<sup>L</sup> = 0.14, β<sup>W</sup> = 0.13. Figure 11c shows the optimal configuration of *C*Q,in is: the Kriging model with a Gaussian kernel and correlation parameters, β<sup>H</sup> = 0.11, β<sup>L</sup> = 0.13, β<sup>W</sup> = 0.17. Figure 11d shows the optimal configuration for the *C*Q,out is: the Kriging model with a Gaussian kernel and correlation parameters, β<sup>H</sup> = 5.52, β<sup>L</sup> = 0.10, β<sup>W</sup> = 2.50. The details are also listed in Table 6. For all hydraulic indices, the performance of the Kriging model is better than that of the other models. The median errors of ξ and *C*<sup>V</sup> are only 0.98% and 0.57% respectively, and the maximum errors are around 5%. Although the maximum error of the *C*Q,out is slightly higher than the others, the value is acceptable as the complex flow separation phenomenon in the outflow condition makes the index more volatile.

**Figure 11.** Surrogate model selection of hydraulic performance indices: (**a**) ξ; (**b**) *C*V; (**c**) *C*Q,in; (**d**) *C*Q, out.

**Table 6.** Optimal surrogate model configurations and their Predictive Estimation of Model Fidelity (PEMF) errors.


#### *5.2. Discussion of Optimized Results of the Inlet*/*Outlet*

Based on the improved process, objective functions are set to find the minimum values for ξ and *C*V. Figure 12 demonstrates the optimization process for finding the optimal shape, in which the red star point indicates the final selection. For the purpose of presenting the optimal result better, the comparisons between the optimal and the original case are listed in Table 7. In the multi-objective optimum situation, ξtotal is reduced by 6.42%, *C*V,total is reduced by 40.28% compared with the original shape.

**Figure 12.** The process of objectives evolved through NSGA-II (dimensions in m). **Table 7.** The comparison of parameters between baseline case and optimum case.


Figure 13 compares velocity contours in the orifices under the inflow condition between the baseline and the optimum design. The inflow velocity contours in the optimum case are found to change little, and the flow velocities in the anti-swirl segment are slightly increased due to the decrease of the optimized inlet height.

**Figure 13.** The comparison of inflow velocity distribution between the baseline and optimum design.

As shown in Figure 14, the change between the baseline and optimum case is more obvious when the water flows out from the tunnel. For the baseline case, the flow separation phenomenon occurs at the upper part of the cross sections in the diffusion segment. Due to insufficient flow diffusion, the separation phenomenon deteriorates into a backflow at the top of the trash rack cross sections, which should be prevented as much as possible in the design. Figure 14a shows the flow separation arises early at *x* = 40 m because of the bad original design. As the flow continues to develop, the flow separation deteriorates and induces reverse velocity at the top of the whole rectification segment and part of the anti-vortex segment.

For the design optimized by NSGA-II, the flow separation region is very close to the walls in the diffusion segment, and the reasonable shape brings in a more desired flow phenomenon. Figure 14b shows the desired flow pattern after optimization. Because of the optimization for the diffusion segment, more uniform velocity distributions are observed in the section plane. The flow separation appears to be slight in the small areas at the junction of the diffusion and rectification segment. Consequently, there is hardly any backflow at the top of the orifices.

**Figure 14.** The comparison of outflow velocity distribution between the baseline and optimum design.

In order to further compare the velocity distribution at the orifice cross-section (*x* = 10), the velocity data from the central line of the section is analyzed. As we can see from Figure 15a, there is little change for the normalized velocity profiles between the baseline and optimum design. Figure 15b demonstrates the velocity profiles along the height direction under the outflow condition. In the baseline case, the main stream of the velocity locates at the lower part, and even negative values appear at the top of the middle orifices. For the optimum design, the velocity distributes more uniformly along the height direction, and the mainstream is in the central section of the height. The maximum value of the normalized velocity is 1.44, which is located at the middle orifice with the corresponding *y*/H value of 0.47. Compared with the original shape, the maximum value is reduced by 49.36%.

**Figure 15.** Velocity profiles at the orifice cross-section: (**a**) inflow; (**b**) outflow.

#### **6. Conclusions**

In this article, the CFD numerical simulation and multi-objective surrogate-based optimization strategy method (COSMOS and NSGA-II) were combined to optimize intake shape effectively. The CFD method applied in this paper was validated by a physical experiment carefully designed to meet bidirectional flow requirements for the inlet/outlet. For the purpose of determining a good compromise between the two modes of operation, reasonable weights were defined to better evaluate the overall performance. Then, suitable surrogate models based on COSMOS were utilized to build functions, and NSGA-II was chosen to complete the optimization.

The optimum shape of the diffusion part in a PHES can be achieved automatically through the whole process, including numerical model building, CFD simulation, optimal surrogate model selection, and multi-objective optimization strategy. There were 60 sampling points generated by the OLHS employed to establish suitable surrogate models based on COSMOS, while the PEMF method was applied to estimate errors. The results indicate that the Kriging with a Gaussian kernel is best for the overall head loss coefficient and the shape parameter β<sup>H</sup> = 2.20, β<sup>L</sup> = 9.41, β<sup>W</sup> = 0.11; the best configuration of the overall velocity distribution coefficient is the Kriging model with an Exponential kernel and β<sup>H</sup> = 1.05, β<sup>L</sup> = 0.14, β<sup>W</sup> = 0.13.

Finally, the NSGA-II algorithm was applied to generate Pareto optimal results, and the final optimal point was selected from the Pareto set through the TOPSIS decision making method. The optimization results show that compared with the original shape the overall head loss coefficient decreases by 6.42% and the overall velocity distribution coefficient decreases by 40.28%. This study demonstrates this new optimization process is a good choice for inlet/outlet engineering designers.

**Author Contributions:** Conceptualization, X.G. and B.S.; methodology, H.Z. (Hongtao Zhu) and Y.T.; software, H.Z. (Han Zhang) and Z.Q.; validation, Z.Q.; formal analysis, H.Z. (Han Zhang); investigation, H.Z. (Hongtao Zhu) and Y.T.; data curation, Z.Q.; writing—original draft preparation, H.Z. (Hongtao Zhu) and Y.T.; writing—review and editing, H.Z. (Hongtao Zhu) and B.S.; visualization, Z.Q.; supervision, X.G. and B.S.

**Funding:** This research was funded by the Science Fund for Creative Research Group of the National Natural Science Foundation of China (51621092), the National Natural Science Foundation of China (51279125) and Tianjin Municipal Natural Science Foundation (15JCYBJC22600).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
