**E**ff**ect of Rotor Spacing and Duct Di**ff**usion Angle on the Aerodynamic Performances of a Counter-Rotating Ducted Fan in Hover Mode**

#### **Woo-Yul Kim, Santhosh Senguttuvan and Sung-Min Kim \***

School of Mechanical Engineering, Sungkyunkwan University, 300 Cheoncheon-dong, Suwon 16419, Korea; samsun215@skku.edu (W.-Y.K.); santhosh@skku.edu (S.S.)

**\*** Correspondence: smkim@skku.edu; Tel.: +82-31-290-7433

Received: 21 August 2020; Accepted: 20 October 2020; Published: 23 October 2020

**Abstract:** The aerodynamic performance of a counter-rotating ducted fan in hover mode is numerically analyzed for different rotor spacings and duct diffusion angles. The design of the counter-rotating fan is inspired by a custom-designed single rotor ducted fan used in a previous study. The numerical model to predict the aerodynamic performance of the counter-rotating ducted fan is developed by adopting the frozen rotor approach for steady-state incompressible flow conditions. The relative angle between the front and the rear rotor is examined due to the usage of the frozen rotor model. The results show that the variation of thrust for the different relative angles is extremely low. The aerodynamic performances are evaluated by comparing the thrust, thrust coefficient, power coefficient, and figure of merit (FOM). The thrust, thrust coefficient, and FOM slightly increase with increasing rotor spacing up to 200 mm, regardless of the duct diffusion angle, and reduce on further increase in the rotor spacing. The duct diffusion angle of 0◦ generates about 9% higher thrust and increases the FOM by 6.7%, compared with the 6◦ duct diffusion angle. The duct diffusion angle is highly effective in improving the thrust and FOM of the counter-rotating ducted fan, rather than the rotor spacing.

**Keywords:** thrust coefficient; power coefficient; figure of merit; frozen rotor; UAV

#### **1. Introduction**

Unmanned aerial vehicles (UAV) are of global interest, since they can perform versatile tasks in the military, search and rescue, agriculture, and transportation fields. Particularly, rotary wing UAVs with vertical take-off and landing (VTOL) propulsion systems are highly desired due to high maneuverability and the ability to take-off and land vertically. In rotary wing UAVs, a ducted fan where a duct surrounds the rotors is more efficient in producing thrust than an open rotor [1]. According to actuator disk theory, the duct reduces the blade tip loss, enabling the ducted fan to double the thrust [2]. The duct also reduces the rotor noise and protects the high-speed rotors from the external environment.

The design of the duct plays a vital role in improving the aerodynamic performance of the ducted fan. One of the key considerations in the design of the ducted fan is the duct profile. The duct profile is usually based on airfoil shapes, due to aerodynamical advantage. Yilmaz et al. [3] experimentally studied the aerodynamic performance of a ducted fan in hover mode for five different duct profiles, based on National Advisory Committee for Aeronautics (NACA) airfoils. Pressure distributions on the inner surface of the duct and velocity profiles at both the duct inlet and outlet planes were investigated for the five different duct profiles. It was found that the propulsive efficiency was significantly affected by the duct profile. Bontempo and Manna [4] numerically investigated the effect of the camber and thickness of the duct profile on the aerodynamic performance of ducted fans. They found that the propulsive efficiency of the duct increased with the increase in the camber and thickness of the

duct profile. Xu et al. [5] developed an asymmetrical duct with a rounded and enlarged leading edge. The leading edge is designed as inflatable so that it can inflate and retract when needed [6,7]. The inflated leading edge is designed to encounter the flow separation that typically occurs around the inner surface of the duct when the ducted fan transits from forward flight (crosswind) to hovering, or from hovering to forward flight. The developed duct model eliminates the flow separation at the crosswind conditions of a 30 m/s velocity at a 50◦ angle of attack. Although the proposed duct model was suitable to deal with the flow separation, it involved certain complications concerned with the practical implementation of the inflatable material in ducted-fan UAVs. Graf et al. [8] conducted experiments to examine the aerodynamic characteristics of a ducted fan for five duct profiles with different leading edge radius of curvatures and duct wall thicknesses, both in hover and forward flight conditions. As a result, the characteristics of the flow separation on the leading edge of the duct in hover and forward flight modes were found to be different. The duct profile that produced the best performance was different for the hover and forward flight modes.

Apart from duct design, a counter-rotating system, where two rotors rotate in the opposite direction on the same axis, can increase the thrust and aerodynamic performance of a ducted fan. The counter-rotating system has been successfully applied to wind turbines [9,10], axial flow pumps [11,12], axial fans [13,14], and propellers [15,16]. A key advantage of this system is that it can produce a higher thrust than a single rotor, and the torque produced by the two rotors is canceled by each other. The essential parameters affecting the aerodynamic performance in the counter-rotating fan are tip clearance and rotor spacing. Ryu et al. [17] examined the effect of the tip clearance of the front and rear rotors on the aerodynamic performances of the counter-rotating ducted fan in hover mode. In the counter-rotating ducted system, the thrust is affected by relative tip clearances of the front and rear rotors. It was found that, in order to improve the aerodynamic performance of the counter-rotating ducted fan, the tip clearance need not be a minimum for both the front and rear rotors. Among various tip clearance configurations in the study, the configuration with a small tip clearance in the front rotor and a large tip clearance in the rear rotor produced the highest thrust coefficient. Moreover, the total thrust of the counter-rotating ducted fan was significantly affected by the tip clearance of the rear rotor. Han et al. [18] conducted experiments and simulations to investigate the aerodynamic performances of a counter-rotating fan with and without ducts for different blade pitch angles, rotor spacings, and tip clearances in hover mode. The results showed that the counter-rotating fan with a duct generated more thrust than the open counter-rotating fan because of the negative pressure region around the duct's leading edge created by the duct. They also found that as the rotor spacing increased, the thrust coefficient and figure of merit of the counter-rotating fan, both with and without the duct, increased. However, when the rotor spacing became higher than the rotor radius, both the thrust coefficient and figure of merit became constant.

Numerical modeling of the counter-rotating fan necessitates the consideration of the relative motion between the two rotors. This relative motion can be modeled by either a mixing plane, a frozen rotor, or a moving mesh approach. Of the three methods, the moving mesh method produces accurate, realistic flow physics by employing the unsteady coupling of the rotors. However, the moving mesh approach requires high computational effort, resulting in researchers and industrial fan designers opting for either the mixing plane or the frozen rotor approach. Besides that, for steady-state simulations, the mixing plane and frozen rotor methods produce reasonably accurate results. A major disadvantage of the mixing plane approach is that it is not capable of predicting the effect of wakes from the top rotor on the downstream rotor of the counter-rotating fan system. Compared to the mixing plane approach, the frozen rotor approach calculates the non-uniform circumferential velocity and pressure distributions to accurately predict the wake mixing in the downstream rotor and rotor-to-rotor flow physics [19].

Despite several studies conducted to enhance the aerodynamic performance of UAVs, there is still scope for improvement. It is highly important to further study and improve the aerodynamic performance of UAVs. In the present study, the effects of rotor spacing and the duct diffusion angle on the aerodynamic performances of a counter-rotating ducted fan in hover mode are investigated numerically using the frozen rotor approach. A ducted fan model, designed by Akturk and Camci [20], is used for the present numerical study due to readily available detailed experimental data. Three essential aerodynamic performance parameters—thrust coefficient, power coefficient, and figure of merit—are examined for different rotor spacings and duct diffusion angles.

#### **2. Numerical Model**

#### *2.1. Model Description*

Figure 1 shows the schematic of the three-dimensional counter-rotating ducted fan based on the fan and duct models of Akturk and Camci [20], where they studied a single ducted fan. The design of both of the fans in the present study was identical to the single fan used by Akturk and Camci [20]. The duct diameter was 559 mm, and the fan had eight blades. The two rotors (fans) had the same tip clearance of 1.71% and rotated in opposite directions. Table 1 summarizes the detailed dimensions of the counter-rotating ducted fan. The rotor spacing (s) was changed from 120 mm to 240 mm in increments of 40 mm. While changing the rotor spacing, the length of duct was also equally changed to maintain the length of the duct diffuser at 117.85 mm. Two duct diffusion angles (θ) of 0◦ and 6◦ were considered. θ

**Figure 1.** Schematic of counter-rotating ducted fan.

**Table 1.** Counter-rotating ducted fan dimensions.


#### *2.2. Governing Equation*

The three-dimensional numerical model to predict the aerodynamic performances of the counter-rotating ducted fan in hover mode was developed using the commercial Computational Fluid Dynamics (CFD) software Ansys CFX 19.1 [21]. The Reynolds-averaged Navier–Stokes (RANS) equations for steady-state, turbulent, and incompressible flow conditions with constant properties were used as the governing equations. The time-averaged continuity and momentum [21] are expressed as follows:

Continuity equation:

$$\frac{\partial \left(\rho u\_j\right)}{\partial \mathbf{x}\_i} = \mathbf{0} \tag{1}$$

Momentum equation for the stationary domain:

$$\frac{\partial \left(\rho u\_{\rangle}\right)}{\partial \mathbf{x}\_{i}} = -\frac{\partial P}{\partial \mathbf{x}\_{i}} + \mu \frac{\partial^{2} u\_{i}}{\partial \mathbf{x}\_{j}^{2}} + \rho g\_{i} \tag{2}$$

Momentum equation for the rotating domain:

$$\frac{\partial \left(\rho u\_{\rangle}\right)}{\partial \mathbf{x}\_{i}} = -\frac{\partial P}{\partial \mathbf{x}\_{i}} + \mu \frac{\partial^{2} u\_{i}}{\partial \mathbf{x}\_{j}^{2}} + \rho g\_{i} + \mathbf{S}\_{M} \tag{3}$$

where the source term, *SM*, includes the centrifugal force and Coriolis force for the rotating reference frame, which is expressed as

$$\mathcal{S}\_M = -\rho \left[ 2\overrightarrow{\boldsymbol{\omega}\_R} \times \overrightarrow{\boldsymbol{u}} + \overrightarrow{\boldsymbol{\omega}\_R} \times \left( \overrightarrow{\boldsymbol{\omega}\_R} \times \overrightarrow{\boldsymbol{r}} \right) \right] \tag{4}$$

In the present study, the shear stress transport (SST) k-ω turbulent model developed by Menter [22,23] was used to solve the turbulent flow field. As per the SST k-ω turbulence model, the boundary layer was calculated by the standard k-ω turbulence model, and the freestream region was calculated by the k-ε turbulence model, which was incorporated by applying a blending function. The SST k-ω turbulence model accurately predicts the adverse pressure gradient and flow separation at the wall [22,23]. The turbulent kinetic energy and dissipation rate are expressed as follows:

Turbulent kinetic energy equation:

$$\frac{\partial \left(\rho u\_{\dot{j}} k\right)}{\partial \mathbf{x}\_{k}} = P - \beta^{\*} \rho \omega k + \frac{\partial}{\partial \mathbf{x}\_{\dot{j}}} \left[ \left( \mu + \sigma\_{k} \mu\_{l} \right) \frac{\partial k}{\partial \mathbf{x}\_{\dot{j}}} \right] \tag{5}$$

Dissipation rate equation:

$$\frac{\partial \left(\rho u\_{\rangle}\omega\right)}{\partial \mathbf{x}\_{\rangle}} = \frac{\gamma}{v\_{\text{f}}} \mathbf{P} - \beta \rho \omega^2 + \frac{\partial}{\partial \mathbf{x}\_{\rangle}} \left[ (\mu + \omega\_{\text{w}} \mu\_{\text{l}}) \frac{\partial k}{\partial \mathbf{x}\_{\rangle}} \right] + 2(1 - \mathbf{F}\_{1}) \frac{\rho \sigma\_{\text{w2}}}{\omega} \frac{\partial k}{\partial \mathbf{x}\_{\rangle}} \frac{\partial \omega}{\partial \mathbf{x}\_{\rangle}} \tag{6}$$

The production limiter (*P)* in Equation (5) is expressed as

$$P = \tau\_{i\bar{j}} \frac{\partial u\_i}{\partial \mathbf{x}\_{\bar{j}}} \tag{7}$$

where

$$
\pi\_{ij} = \mu\_l \left( 2S\_{ij} - \frac{2}{3} \frac{\partial u\_k}{\partial \mathbf{x}\_k} \delta\_{ij} \right) - \frac{2}{3} \rho k \delta\_{ij} \tag{8}
$$

$$S\_{ij} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \tag{9}$$

and the turbulent eddy viscosity is expressed as

$$\mu\_{\rm f} = \frac{\rho a\_1 k}{\max(a\_1 \omega, \,\Omega F\_2)} \tag{10}$$

Let φ<sup>1</sup> represent the constants σ*k*<sup>1</sup> , σω1, and β<sup>1</sup> in the k-ω turbulence model, φ<sup>2</sup> represent the constants σ*k*<sup>2</sup> , σω2, and β<sup>2</sup> in k-ε turbulence model, and φ represent the constants σ*<sup>k</sup>* , σω, and β in the k-ω SST turbulence model, respectively. The relationship between φ, ]φ1, and φ<sup>2</sup> is expressed as

$$
\phi = F\_1 \phi\_1 + (1 - F\_1)\phi\_2 \tag{11}
$$

The additional functions in the SST k-ω turbulence model are as follows:

$$F\_1 = \tanh(\arg\_1^4) \tag{12}$$

$$\arg\_1 = \min \left[ \max \left( \frac{\sqrt{k}}{\beta^\* \omega d}, \frac{500v}{d^2 \omega} \right), \frac{4\rho \sigma\_{\omega 2} k}{\mathbb{C}D\_{k\omega} d^2} \right] \tag{13}$$

$$\text{CD}\_{k\omega} = \max \left( 2\rho \sigma\_{a2} \frac{1}{\omega} \frac{\partial k}{\partial \mathbf{x}\_{\parallel}} \frac{\partial \omega}{\partial \mathbf{x}\_{\parallel}}, 10^{-20} \right) \tag{14}$$

$$F\_2 = \tanh(\arg\_2^2) \tag{15}$$

$$\arg \mathbf{g}\_2 = \max \left( 2 \frac{\sqrt{k}}{\beta^\* a \alpha d}, \frac{500v}{d^2 \omega} \right) \tag{16}$$

where *d* is the distance from the field point to the nearest wall and Ω = p 2*WijWij*, with

$$\mathcal{W}\_{ij} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} - \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \tag{17}$$

All the constants are as follows:

$$\begin{aligned} \gamma\_1 = \frac{\beta\_1}{\beta^\*} - \frac{\sigma\_{\omega 1} \mathbf{x}^2}{\sqrt{\mathfrak{p}^\*}}, \; \gamma\_2 = \frac{\beta\_2}{\beta^\*} - \frac{\sigma\_{\omega 2} \mathbf{x}^2}{\sqrt{\mathfrak{p}^\*}}, \; \sigma\_{k1} = 0.85, \; \sigma\_{\omega 1} = 0.5, \; \beta\_1 = 0.075, \; \sigma\_{k2} = 1.0, \\\sigma\_{\omega 2} = 0.856, \; \beta\_2 = 0.0828, \; \beta^\* = 0.09, \; \kappa = 0.41, \; a\_1 = 0.31. \end{aligned} \tag{18}$$

#### *2.3. Numerical Method*

The numerical simulations were performed as three-dimensional with steady-state and incompressible flow assumptions. Figure 2 shows the grid system, boundary, and interface conditions of the counter-rotating ducted fan. The computational domain consisted of a stationary domain, surrounding the duct and part of the hub, and two rotating domains surrounding each rotor and the remaining part of the hub. The unstructured grid was used for both the stationary and the rotating domain, and the structured hexahedral grid was used for the external fluid domain. The dimensionless distance from the wall to the first node of the mesh, defined as y+, is crucial in the SST k-ω turbulence model. The inflation layers were used to achieve a y+ value less than 2.5 along all the walls to capture fluid interactions in the viscous sublayer. The axial length of the computational domain was 3.0 D and 6.0 D from the origin, located at the center of the front rotor, and the radial length was 2.5 D from the *y*-axis. The multiple reference frame (MRF) method was used to simulate the rotation of the rotors. The rotating domain was considered to rotate relative to the stationary domain by adding the centrifugal force and Coriolis force (see Equations (2) and (3)). The frozen rotor, mixing plane (stage), and sliding mesh methods were the available methods to model the interface between rotating and stationary domains in the MRF method. In the present study, the frozen rotor method was used for the interface between the rotating and stationary domains, as well as the interface between the two rotating domains. When using the frozen rotor, the rotor in the rotating domain is considered to be in a frozen state; this means the rotor position is fixed relative to the stationary domain. In the frozen rotor method, the results may vary depending on the rotor position. Hence, in the present study, the positions of the front and rear rotors were tested for four relative angles: 0◦ , 11.25◦ , 22.5◦ , and 33.75◦ . The opening boundary condition with zero-gauge pressure was applied to the external surfaces of the fluid domain to allow inward or outward airflow in hover mode. The no-slip condition was applied to the walls of the duct, hub, and rotors, except for the inner wall, the duct adjacent to the rotating domains. The counter-rotating wall boundary condition with an equal velocity in the opposite direction corresponding to the rotating domain was applied to the inner wall of the duct. To reduce the computational cost, only one-eighth of the counter-rotating ducted fan was modeled in the computational domain by applying the periodic boundary conditions to the side surfaces.

**Figure 2.** The grid system, boundary, and interface conditions of the counter-rotating ducted fan.

#### *2.4. Grid Independence Test and Validation*

The aerodynamic performance parameters, thrust coefficient, power coefficient, and figure of merit are expressed as follows:

24

35

Thrust coefficient:

$$C\_T = \frac{\text{Thrust}}{\rho N^2 D^4} \tag{19}$$

Power coefficient:

$$C\_P = \frac{\text{Power}}{\rho N^3 D^5} \tag{20}$$

Figure of merit (FOM):

$$FOM = \frac{\mathbf{C}\_T^{3/2}}{\sqrt{2}\mathbf{C}\_P} \tag{21}$$

√2 The single ducted fan model of Akturk and Camci [20] was simulated to conduct grid independence testing and validation. The grid independence test was conducted using four different grid systems comprising 3.98, 4.79, 5.87, and 7.02 million cells at 1500 rpm. The thrust coefficient of the different grid systems is compared in Figure 3a. The discrepancy in the thrust coefficient among the 4.79, 5.87, and 7.02 million cell grids was less than 1.0%, and, thus, the grid system with 4.79 million cells was selected. The total simulation time for the selected mesh size was about 10 h in a machine with an Intel Core i7 6700 (4 physical cores and 8 threads). For validation, the single ducted fan was simulated under six different rotational speeds, corresponding to the experiments of Akturk and Camci [20]. Figure 3b compares the thrust generated for the six different rotational speeds in the present numerical study with the experiment and numerical results of Akturk and Camci [20]. The present numerical results show good agreement with the results of Akturk and Camci [20].

**Figure 3.** (**a**) The thrust coefficient for different grid systems, and (**b**) comparison of thrust on single rotor ducted fan of the present numerical study with the experiment and numerical results of Akturk and Camci [20].

#### **3. Results and Discussion**

#### *3.1. E*ff*ect of Relative Angle between Front and Rear Rotor*

Initially, the effect of the relative angle between the front and rear rotor was tested for four relative angles: 0◦ , 11.25◦ , 22.5◦ , and 33.75◦ . Each relative angle configuration was simulated for three rotational speeds—1500 rpm, 2100 rpm, and 2700 rpm—and two rotor spacings of 120 mm, and 200 mm. The diffusion angle of the duct for all of the relative angle configurations was 6◦ . Figure 4 shows the thrust generated at each relative angle configuration for the different rotational speeds and rotor spacings. The thrust was not significantly affected by the relative angle at 1500 rpm for both rotor spacings. However, minor variations in the thrust were observed as the rotational speed increased for both rotor spacings.

**Figure 4.** Thrust generated at four relative angles between the front and rear rotors for different rotational speeds and rotor spacings.

Figure 5a–d shows the total pressure contours at Plane 3 for the different relative angles under the conditions of a 2700 rpm rotational speed and 200 mm rotor spacing. Refer to Figure 6 for the location of Plane 3 in the domain. The pressure contours of the 2700 rpm rotational speed and 200 mm rotor spacing cases were compared since the thrust was slightly affected by the relative angles as the rotational speed increased. The high-pressure region, widely distributed across the center of the interface, shifted periodically with the change in relative angle. The difference in the pressure distribution for the different relative angles was negligible. Hence, the effect of the relative angle between the front and rear rotor was negligible, and the 0◦ relative angle was considered for further simulations to evaluate the aerodynamic performance of the counter-rotating ducted fan. Similar relative angle studies conducted in a gas turbine pre-swirl system [24], an axial turbine [25], and a centrifugal pump [26] also concluded that the effect of the relative angle is negligible.

θ ω **Figure 5.** Total pressure contour and four relative angles between the front and rear rotors of (**a**) 0◦ ; (**b**) 11.25◦ ; (**c**) 22.5◦ ; and (**d**) 33.75◦ at Plane 3, when s = 200 mm, θ = 6 ◦ , and ωR = 2700 rpm. Refer to Figure 6 for the location of Plane 3. θ ω

(1) Plane 1 - Interface plane between front rotating domain and stationary domain (1) Plane 1 - Interface plane between front rotating domain and stationary domain

(2) Plane 2 - Interface plane between front rotating domain and rear rotating domain (2) Plane 2 - Interface plane between front rotating domain and rear rotating domain

(3) Plane 3 - Interface plane between rear rotating domain and stationary domain (3) Plane 3 - Interface plane between rear rotating domain and stationary domain


**Figure 6.** Plane locations used for plotting.

#### *3.2. Aerodynamic Performance Analysis of Counter-Rotating Ducted Fan*

Numerical simulations were done for different rotor spacings and duct diffusion angles (see Section 2.1) under five rotating speeds of 1500 rpm, 1800 rpm, 2100 rpm, 2400 rpm, and 2700 rpm in hover mode. Figure 7 shows the thrust for the different rotor spacings and duct diffusion angles. The cases with a 0◦ duct diffusion angle generated about 9% higher thrust than the cases with a 6 ◦ duct diffusion angle. The rotor spacing variation had a minimal effect on improving the thrust. The maximum thrust for both duct diffusion angles was observed for the rotor spacing of 200 mm, which can also be verified quantitatively from Table 2, where the thrust, thrust coefficient, and the power coefficient are shown for all cases. The thrust increased with increasing rotor spacing up to 200 mm, and then it started to reduce upon further increase of the rotor spacing. The most substantial thrust increment occurred when the rotor spacing increased from 120 mm to 160 mm. The thrust of the 200 mm rotor spacing increased by about 1.3–1.5% compared with that of the 120 mm rotor spacing. Like the tendency of the thrust results, the maximum thrust coefficient was observed when the rotor spacing was 200 mm. Moreover, the thrust coefficient was higher when the duct diffusion angle was 0◦ rather than 6◦ . However, the power coefficient was at a minimum when the rotor spacing was 120 mm and, in terms of the diffusion angle, the minimum power coefficient was attained for a 6◦ diffusion angle. In summary, from Table 2, it can be observed that the thrust, thrust coefficient, and power coefficient became higher either when the duct diffusion angle decreased from 6◦ to 0◦ or when the rotor spacing increased to 200 mm. Figure 8 shows the figure of merit for different rotor spacings and duct diffusion angles. The trend of FOM is similar to the thrust results. The FOM increased by about 1% with increasing rotor spacing up to 200 mm from 120 mm, and it reduced gradually upon a further increase of rotor spacing. The duct diffusion angle was more effective than the rotor spacing for improving the FOM. The FOM increased by about 6.7% when the duct diffusion angle decreased from 6 ◦ to 0◦ .


**Table 2.** Thrust, thrust coefficient, and power coefficient of all the counter-rotating ducted fan cases.

**Figure 7.** Thrust for different rotor spacings and duct diffusion angles.

**Figure 8.** Figure of merit for different rotor spacings and duct diffusion angles.

θ θ Since the rotor spacing had a minimal effect on improving the aerodynamic performance of the counter-rotating ducted fan, henceforth, only the results of the cases with rotor spacings of 120 mm and 200 mm for θ = 0 ◦ , 6◦ simulated at the maximum rotational speed of 2700 rpm are discussed. Figure 9a–d displays the total pressure contour in Plane 3. See Figure 6 for the location of Plane 3.

The negative pressure region close to the duct shows rotor tip leakage loss. As the rotor spacing increases, the thrust increases by about 1%, and, hence, the positive pressure area is extended marginally. The change in the positive pressure area as the rotor spacing increases is not significant. However, when the duct diffusion angle reduces from 6◦ to 0◦ , the magnitude of the positive pressure increases significantly.

Figure 10a–d shows the axial velocity contours in Plane 4. See Figure 6 for the location of Plane 4. The axial velocity near the hub increases in the radial direction due to the 9% thrust increase as the duct diffusion angle reduces from 6◦ to 0◦ . As the duct diffusion angle is increased to 6◦ , the cross-sectional area of the duct is also increased, resulting in a lower axial velocity. The high axial velocity region is widened for an increase in the rotor spacing and a decrease in the duct diffusion angle.

 s = 120 mm Figure 11a–d shows the axial velocity profiles in the radial direction at different planes. See Figure 6 for the location of the planes. The radial length is normalized with the duct radius. In Figure 11a, the magnitude and profile of the axial velocities located at Plane 1 for all cases are fairly similar. Figure 11b indicates that the axial velocity increases at Plane 2 when the rotor spacing is narrow, which is caused by the influence of the rear rotor. The axial velocities in Plane 3 show variation only after the normalized radial length of about 0.7 (see Figure 11c). Moreover, at Plane 3, the magnitude of the axial velocity is higher when s = 200 mm, due to high thrust. In Figure 11d, the magnitude of the axial velocity at Plane 5 for θ = 6 ◦ is reduced, resulting from the increased cross-sectional area of the duct.

s = 120 mm

20 60 80

0.55 0.56

0.58

0.57

0.59

0.60

0.61

**Figure of merit [-]**

0.63

0.62

0.65

0.64

120

100

**Thrust [N]**

160

140

200

180

s = 200 mm s = 160 mm s = 240 mm

θ = 0°

θ = 6° s = 120 mm s = 160 mm s = 200 mm s = 240 mm

**Rotational speed [rpm]** 1200 1500 1800 2100 2400 2700 3000

> θ = 6° s = 120 mm s = 160 mm s = 200 mm s = 240 mm

θ = 0°

**Rotational speed [rpm]** 1200 1500 1800 2100 2400 2700 3000

s = 120 mm

s = 200 mm s = 160 mm s = 240 mm s = 120 mm

s = 200 mm s = 160 mm s = 240 mm

 **Figure 9.** Total pressure contour in Plane 3 at 2700 rpm for: (**a**) θ = 6 ◦ , s = 120 mm; (**b**) θ = 6 ◦ , s = 200 mm; (**c**) θ = 0 ◦ , s = 120 mm; and (**d**) θ = 0 ◦ , s = 200 mm. See Figure 6 for the location of Plane 3. θ

 **Figure 10.** Axial velocity contour in Plane 4 at 2700 rpm for (**a**) θ = 6 ◦ , s = 120 mm; (**b**) θ = 6 ◦ , s = 200 mm; (**c**) θ = 0 ◦ , s = 120 mm; and (**d**) θ = 0 ◦ , s = 200 mm. See Figure 6 for the location of Plane 4.

**Figure 11.** Axial velocity profile along the radial direction at 2700 rpm in (**a**) Plane 1; (**b**) Plane 2; (**c**) Plane 3; and (**d**) Plane 5. See Figure 6 for the location of the planes.

#### **4. Conclusions**

The present study explores the aerodynamic performances of a counter-rotating ducted fan in hover mode numerically. The effect of different rotor spacings and duct diffusion angles on the aerodynamic performances is examined. Key findings from this study are as follows:


### **5. Future Work**

The present study discussed the effects of rotor spacing and the duct diffusion angle on the aerodynamic performances of a counter-rotating ducted fan in hover mode. However, to enhance the aerodynamic performance of the UAV, there are still several other design parameters to be considered. Different duct lip shapes and tip clearances could be considered to optimize the counter-rotating ducted fan toward the development of high-performance, efficient UAVs. Moreover, apart from the frozen rotor approach, it is critical to explore various other numerical approaches, such as harmonic analysis techniques, to numerically model the counter-rotating ducted fan.

**Author Contributions:** Validation, methodology, software, writing-original draft preparation, W.-Y.K.; validation, writing—original draft preparation, S.S.; conceptualization, methodology, supervision, writing—review and editing, S.-M.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 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/).

### *Article* **Thermal Performance of T-shaped Obstacles in a Solar Air Heater**

### **Seung-Yong Ahn <sup>1</sup> and Kwang-Yong Kim 2,\***


Received: 18 September 2020; Accepted: 13 October 2020; Published: 17 October 2020

**Abstract:** This paper proposes T-shaped ribs as obstacles attached to the heat absorber plate in a rectangular solar air heater to promote heat transfer. The thermal and aerodynamic performance of the solar heater was numerically evaluated using three-dimensional Reynolds-averaged Navier–Stokes equations with the shear stress transport turbulence model. A parameter study was performed using the ratios of rib height to channel height, rib width to channel width, and rib width to rib height. The area-averaged Nusselt number and friction factor were selected as the performance parameters of the solar air heater to evaluate the heat transfer and friction loss, respectively. In addition, the performance factor was defined as the ratio of the area-averaged Nusselt number to the friction factor. The maximum area-averaged Nusselt number was found at h/e = 0.83 for a fixed rib area. Compared with triangular ribs, the T-shaped ribs showed up to a 65 % higher area-averaged Nusselt number and up to a 49.7% higher performance factor.

**Keywords:** solar air heater; ribs; Nusselt number; friction factor; Reynolds-averaged Navier–Stokes equations

#### **1. Introduction**

Recently, as the problem of environmental pollution due to the use of fossil fuels has emerged, interest in solar heat systems as a renewable energy source is increasing. A solar air heater (SAH) is a device that heats air flowing over an absorber plate by using solar energy. It is used for space heating, drying of agricultural products, and dehydration of industrial products [1].

Compared to liquid solar heaters, SAHs have a disadvantage of low heat transfer rate, because the density of air is lower than that of a liquid. The heat transfer performance of an SAH is determined by various factors such as the velocity of the flow, the length and depth of the heater, and the shape of the heat absorber plate. The ratio between the area of the actual heat absorber and the absorber area normal to the solar radiation, which is called the absorber shape factor, is an important parameter in the design of SAHs.

Obstacles attached to the heat absorber plate are generally known to increase the heat transfer by enlarging the heat transfer area and enhancing turbulence intensity. However, as the area of the obstacles attached to the heat absorber plate increases, the pressure drop through the SAH also increases. Therefore, many studies have been conducted on the shape of the obstacles to enhance the overall heat transfer performance by increasing the heat transfer while minimizing the pressure drop [2–14].

Kabeel and Mecarik [2] studied the effect of the shape of the heat absorber plate on the performance of an SAH. It was confirmed that the heat transfer performance of the SAH with triangular obstacles was higher than that with longitudinal pins. Moummi et al. [3] performed an energy analysis of an SAH with rectangular plate fins. The results for the heat transfer coefficient were compared with the results obtained for the SAH without obstacles. The collector efficiency factor of the SAH with plate pins increased by 30% compared to the case without obstacles. Essen [4] performed energy and exergy analysis through experiments on obstacles of various shapes attached to an SAH. He concluded that the heat and exergy efficiencies of the SAH with obstacles increased as compared to the SAH without obstacles and the collector efficiency factor depended on the shape, area, orientation, and arrangements of the obstacles. Saini and Saini [5] conducted an experimental study on the effect of the rib shape on the heat transfer performance of an SAH with arc-shaped ribs. They developed correlations for the Nusselt number and friction factor in terms of Reynolds number (Re), relative roughness height, and arc angle of ribs.

Ozgen et al. [6] conducted an efficiency evaluation of an SAH with the flows flowing over the upper and lower surfaces of a heat absorber plate where aluminum cans were attached. Depaiwa et al. [7] experimentally studied the forced convective heat transfer and friction loss for the turbulent flow in an SAH with rectangular winglet vortex generators. They tested an SAH with 20 winglet vortex generators at Reynolds numbers ranging from 5000 to 23,000. The heat transfer rates in this SAH were 174% to 182% higher than that of the SAH with a smooth absorber plate. Bekele et al. [8] performed an experimental and numerical analysis on the effects of the height and longitudinal pitch of triangular obstacles on the heat transfer in an SAH. Through the analysis of the internal flow of the SAH, it was reported that recirculating flows occurring around the triangular obstacles effectively increased the turbulence intensity, thereby increasing the heat transfer rate by 3.6 times compared to the SAH without obstacles. Yadav et al. [9] conducted an experimental study on the heat transfer performance and friction factor of the turbulent flow over a heat absorber plate with circular protrusions arranged in angular arc. The Nusselt number and friction factor were found to be 2.89 times and 2.93 times higher than those of the unobstructed SAH, respectively.

Kulkarni and Kim [10] performed a numerical analysis to find the optimal shape of obstacles attached to an SAH. Numerical analysis was performed for four different shapes and three different arrangements of the obstacles. The Nusselt number and friction coefficient were greatly influenced by the shape and arrangement of the obstacles, and pentagonal obstacles showed the highest performance factor among the tested shapes. Alam and Kim [11] performed a numerical analysis to predict the heat transfer performance of an SAH with conical protrusion ribs. The maximum thermal efficiency of the collector was reported to be 69.8%. They developed correlations for the Nusselt number and friction factor in terms of Reynolds number, relative roughness height, and relative rib pitch. Compared with the results of numerical analysis, the correlations for the Nusselt numbers and friction factor showed average absolute standard deviations of 2.78% and 5.25%, respectively. In addition, studies on SAHs with various types of obstacles such as V-shaped ribs, diamond-shaped ribs, and arc-shaped ribs have been also conducted to date [12–14].

In this study, T-shaped ribs are newly proposed to further improve the performance of an SAH by considering both the heat transfer and pressure drop. The heat transfer and pressure drop characteristics of the SAH with ribs were analyzed using three-dimensional Reynolds-averaged Navier–Stokes (RANS) equations. A parametric study was conducted to confirm the effects of geometric parameters of the T-shaped ribs on the heat transfer and pressure drop in the SAH. The performance of the SAH with the proposed ribs was compared to SAHs with previously developed obstacles.

#### **2. Solar Air Heater Model**

Figure 1 shows the computational domain of the SAH and geometric parameters of the T-shaped rib. Only a half of the entire SAH domain is included in the computational domain using symmetric conditions. The computational domain consists of three sections; inlet Section (800 mm × 150 mm × 50 mm (L1 × W/2 × H)), test Section (1200 mm × 150 mm × 50 mm (L2 × W/2 × H)), and outlet Section (500 mm × 150 mm × 50 mm (L3 × W/2 × H)). This computational domain for evaluating the performance of the SAH was set according to the guidelines of the ASHRAE standard 93–97 [15], and the lengths of the inlet and outlet sections of the domain are 5(WH)0.5 and 2.5(WH)0.5, respectively.

The upper wall (i.e., heat absorber plate) of the test section absorbs solar energy, and air is heated as it flows from the inlet to the outlet of the channel. The obstacles, T-shaped ribs, are attached to the heat absorber plate. Thirteen rows of ribs are arranged in a zigzag, as shown in Figure 1a. The thickness of the rib is 0.5 mm, and t is 8 mm. Each row contains one or one and a half ribs in the computational domain. The first row of ribs is located 71.75 mm downstream of the inlet of the test section.

**Figure 1.** Computational domain of the solar air heater with T-shaped ribs. (**a**) Computational domain, (**b**) Test section.

#### **3. Numerical Methods**

ε ω ω ε In this study, ANSYS-CFX 15.0 [16], a commercial code which employs an unstructured grid, was used to analyze the flow and heat transfer using RANS equations. The shear stress transport (SST) [17] model was used for the analysis of turbulence. The SST model was designed to take advantage of the k-ε and k-ω models by combining these two models using a weighting function so that the k-ω model is applied near the wall and the k-ε model is applied in the other area. The SST model is known to be effective in predicting the flow recirculation due to separation.

Air at 25 ◦C was used as the working fluid, and velocity and static pressure conditions were assigned to the inlet and outlet boundaries, respectively. A constant heat flux condition of 815 W/m<sup>2</sup> was applied to the surface of the heat absorber plate, and the other solid surfaces were assumed to be adiabatic. The symmetric conditions were used at the symmetric plane of the computational domain, and a no-slip condition was used at the wall boundaries. An unstructured tetrahedral grid system

was constructed in most of the computational domain, but prism meshes were placed near the wall to resolve the laminar sublayer. In order to use low-Re modeling for near-wall turbulence, the y+ values of the first grid points from the wall were kept at less than 1.0.

To find the grid dependency of the numerical solution, a test was performed according to the procedure proposed by Roache [18] and Celik and Karatekin [19]. This test analyzes the grid convergence index (GCI), which represents numerical uncertainty through the estimation of discretization error based on Richardson extrapolation. Table 1 shows the results of the GCI analysis where the grid refinement factor (r) was set to 1.3 for three different grid systems (N1, N2, and N3). This test was performed on the SAH with the reference ribs described in Table 2. As the number of grid nodes increases, the Nusselt number tends to converge gradually. When N<sup>2</sup> is used, the extrapolated relative error (*e* 21 *ext*) is 0.0015% and the relative error is 0.0019%, which confirms a relatively small numerical uncertainty. Therefore, based on this result, the optimum grid system is selected as N2, which is shown in Figure 2. ௫௧ ଶଵ

**Table 1.** Results of grid convergence index (GCI) analysis using Richardson extrapolation.


**Table 2.** Ranges and reference values of geometric parameters.


**Figure 2.** Example of grid system.

In order to determine the convergence of the numerical solution, the root mean square residual was reduced to 10−<sup>6</sup> or less, and about 1200 iterations were performed. The computational time required for a single analysis was about 14 h when a personal computer with an Intel Core i7 3.41 GHz CPU was used for the computations.

#### **4. Geometric and Performance Parameters**

In order to examine the effects of the rib shape shown in Figure 1 on the heat transfer performance and pressure drop, three dimensionless parameters were selected; the ratios of the rib height to the channel height (h/H), the rib width to the channel width (e/W), and the rib height to the rib width (h/e). The ranges and reference values of these parameters are shown in Table 2.

Three performance parameters are defined to evaluate the heat transfer and pressure drop in the SAH. The performance parameter related to the heat transfer is defined as follows:

$$F\_{\rm Nu} = \frac{\rm Nu\_o}{\rm Nu\_s} \tag{1}$$

*Nu<sup>o</sup>* is the area-averaged Nusselt number on the heat transfer surface.

$$Nu\_o = \frac{q\_o D\_h}{k\_a (T\_p - T\_s)} \tag{2}$$

where *T<sup>p</sup>* is the average temperature at the surface of the heat absorber plate including the obstacles, *T<sup>s</sup>* is the bulk temperature of the working fluid, *k<sup>a</sup>* is the thermal conductivity of the fluid, *q<sup>o</sup>* is the heat flux at the heat absorber plate, *D<sup>h</sup>* is the hydraulic diameter of the flow channel, and *Nu<sup>s</sup>* is the Nusselt number for a fully developed flow in a smooth channel without obstacles, which is calculated from the following Dittus–Boelter equation.

$$Nu\_{\rm s} = 0.024 Re^{0.8} Pr^{0.4} \tag{3}$$

where Reynolds number, *Re*, is defined using a hydraulic diameter, and *Pr* indicates the Prandtl number.

The definition of the performance parameter related to the pressure drop is as follows:

$$F\_f = \left(\frac{f\_o}{f\_s}\right)^{1/3} \tag{4}$$

Here, *f<sup>o</sup>* is the friction factor in the flow channel with obstacles.

$$f\_o = \frac{2(\Delta P)D\_h}{4\rho L \mathcal{U}^2} \tag{5}$$

where ∆*p* is the pressure drop in the test section, ρ is the density of the working fluid, *U* is the average velocity of the flow, *L* is the length of the flow channel (test section), and *f<sup>s</sup>* is the friction factor for the fully developed flow in the smooth channel, which is calculated from the following modified Blasius equation.

$$f\_s = 0.085 \text{Re}^{-0.25} \tag{6}$$

The performance factor, *PF* [20], for evaluating the performance of an SAH considering both the heat transfer performance and pressure drop at the same time, is defined as follows:

$$PF = \frac{F\_{Nu}}{F\_f} \tag{7}$$

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

To prove the validity of the numerical analysis, the numerical results for the Nusselt number, Nuo, and friction factor, *f* <sup>o</sup>, are compared with corresponding empirical formulas for different Reynolds numbers in the flow channel without obstacles, as shown in Figure 3. The numerical results show good overall agreement with the empirical formulas. As Reynolds number increases, the relative error tends to decrease in the case of Nuo. Compared with the experimental data, the average relative error of the computed Nusselt numbers within the tested Reynolds number range is less than 3.29%, and it is less than 2.71% for the friction factor.

**Figure 3.** Validation of numerical results for smooth duct.

In addition, the validation test was also performed for the case with obstacles. The numerical results are compared with the experimental data obtained by Bekele et al. [8] for the SAH with triangular obstacles under the same boundary conditions, as shown in Figure 4. The numerical results agree well qualitatively with the experimental data. In the results for the Nusselt number, the computational and empirical slopes according to the Reynolds number are exactly the same but, quantitatively, they show a small difference unlike the case without obstacles. The average relative error for the Nusselt number compared to the experimental data within the tested Reynolds number range is less than 4.43%, and that for the friction factor is less than 1.94%, which indicates that the accuracy of the numerical analysis results is acceptable for further calculations.

**Figure 4.** Validation of numerical results for a solar air heater (SAH) with triangular ribs using experimental data of Bekele et al. [8].

≤ ≤

≤ ≤

The three dimensionless geometric parameters of the SAH with proposed T-shaped ribs presented in Table 2 were used to find the effects of the rib shape on the performance parameters, i.e., FNu, F*<sup>f</sup>* , and PF. The parametric study was performed at a Reynolds number of 22,600. In the case of two variables, h/H and e/W, the reference value in Table 2 was applied to the parameter that was not changed, and the area of the rib was changed accordingly during the parametric study. However, in the case of h/e, the test was conducted with the rib area fixed at 375 mm<sup>2</sup> . Therefore, in this case, both the absolute values of h and e changed to keep the area constant according to the change of h/e. In all cases, the size of t was kept constant.

The results of the parametric study for the ratio of the rib height to the channel height (h/H) are shown in Figures 5–9. Figures 5 and 6 show the changes in the performance parameters, FNu and F*f* , defined by Equations (1) and (2), respectively, in a range of 0.3 ≤ h/H ≤ 0.7. The computational values are presented at five points indicated by circular symbols, and the line is a curve fit of the values at these points. It can be seen that FNu has the maximum value between h/H = 0.5 and 0.6 and it gradually decreases after that (Figure 5). On the other hand, F*<sup>f</sup>* continues to increase as h/H increases, as shown in Figure 6. Therefore, as the rib height increases with the fixed rib width, there exists the maximum value in the heat transfer rate, but the pressure drop increases continuously. Similar heat transfer performances are found at h/H = 0.5 and 0.6 (Figure 5), but the pressure drop is lower at the lower rib height (h/H = 0.5), as shown in Figure 6. Therefore, in the PF distribution shown in Figure 7, the maximum performance factor appears at h/H = 0.5.

**Figure 5.** Variation of the performance parameter related to heat transfer (FNu) with the ratio of the rib height to the channel height (h/H) (the ratio of the rib width to the channel width (e/W) = 0.1).

ƒ **Figure 6.** Variation of the performance parameter related to the pressure drop (F*<sup>f</sup>* ) with h/H (e/W = 0.1).

ƒ

**Figure 7.** Variation of performance factor (PF) with h/H (e/W = 0.1).

**Figure 8.** Streamlines on x-y plane (z = 0) between fifth and seventh ribs (dark region in Figure 1b) for different h/H. (**a**) h/H = 0.3, (**b**) h/H = 0.5, (**c**) h/H = 0.7.

**Figure 9.** Nusselt number distributions on the heated surface for different h/H. (**a**) h/H = 0.3, (**b**) h/H = 0.5.

Figure 8 shows the streamlines between the 5th and 7th ribs (indicated by a dark region in Figure 1b) at the symmetric plane (i.e., x-y plane at z = 0) for different h/H. Due to the geometric shape of the rib, two recirculation regions are formed in the direction of the rib height. As h/H increases, the size of the recirculation region near the upper wall increases. Additionally, as h/H increases from 0.3 to 0.5, the reattachment distance rapidly decreases. This rapid reduction in the reattachment distance in this h/H range greatly enhances the heat transfer, as shown in Figure 5, due to early re-development of the thermal boundary layer. In Figure 9, the local Nusselt number distributions on the heat transfer surface for h/H = 0.3 and 0.5 are different. In the case of h/H = 0.5, the Nusselt number level appears to be higher around and downstream of each rib compared to the case of h/H = 0.3.

Figures 10–14 show the effects of the ratio of the rib width to the channel width (e/W) on the performance of the SAH. Figures 10 and 11 show the variations of the two performance parameters, FNu and F*<sup>f</sup>* , respectively, with e/W, in a range of 0.067 ≤ e/W ≤ 0.133. As e/W increases, both FNu and F*<sup>f</sup>* increase almost linearly. This is presumed to be a phenomenon that occurs because both the turbulence intensity and pressure loss increase due to the increase in the area of the obstacles as e increases while the height h is kept constant. In the PF distribution shown in Figure 12, the maximum value is found around e/W = 0.08, but the overall variation with e/W is very small compared to that in Figure 7 for h/H. ≤ ≤ ≤ ≤

Figure 13 shows the streamlines between the 5th and 7th ribs at the symmetric plane for different h/H. It can be seen that as e/W increases, the recirculation region located near the top of the rib increases. Since the reattachment distance does not change with e/W, the heat transfer enhancement according to the change in e/W does not seem to be related to the reattachment distance. Figure 14 shows the distributions of the Nusselt number on the heat transfer surface between the 5th and 7th ribs. As e/W increases, the increased width of the rib increases the width of the high Nusselt number area in the lateral direction, thereby increasing the overall Nusselt number.

**Figure 10.** Variation of FNu with e/W (h/H = 0.5).

**Figure 11.** Variation of F*f*ƒwith e/W (h/H = 0.5).

ƒ

**Figure 12.** Variation of PF with e/W (h/H = 0.5).

(**a**)

(**c**)

**Figure 13.** Streamlines on x-y plane (z = 0) between fifth and seventh ribs for different e/W. (**a**) e/W = 0.067, (**b**) e/W = 0.10, (**c**) e/W = 0.133.

**Figure 14.** Nusselt distributions on the heated surface between fifth and seventh ribs for different e/W. (**a**) e/W = 0.067; (**b**) e/W = 0.10; (**c**) e/W = 0.133.

The variation of FNu with the ratio of the rib height to the rib width (h/e) keeping the rib area constant is shown in Figure 15. At h/e = 0.83, the heat transfer is maximized. This variation of FNu is similar to the variation with h/H shown in Figure 5. However, unlike Figure 5, the variation curve of FNu shown in Figure 15 is more symmetrical as the area of the rib is kept constant in this case. The variation of F*<sup>f</sup>* with h/e shown in Figure 16 is very small compared to the variations with h/H and e/W shown in Figures 6 and 11, respectively. In the cases with h/H and e/W, the area of the rib changes, but in this case with h/e, the area of the rib remains constant. Accordingly, it can be seen that the pressure drop is greatly influenced by the area of the rib rather than the aspect ratio of the rib.

**Figure 15.** Variation of FNu with the ratio of the rib height to the rib width (h/e) (rib area = 375 mm<sup>2</sup> ).

ƒ **Figure 16.** Variation of F*<sup>f</sup>* with h/e (rib area = 375 mm<sup>2</sup> ).

ƒ Since the change in the pressure drop appears small over the entire range of h/e, the variation of the performance factor PF with h/e, shown in Figure 17, shows a qualitatively similar variation to that of FNu (Figure 15). The maximum value of PF also appears at h/e = 0.83.

Figure 18 shows the Nusselt number distributions on the heat transfer surface between the 5th and 7th rib rows. At h/e = 0.83, where the highest heat transfer occurs, the area of the high Nusselt number region is the largest. This indicates that the heat transfer enhancement due to the production of turbulent kinetic energy becomes most effective when h and e have similar sizes while keeping the rib area constant.

To demonstrate the superiority of the heat transfer performance of the SAH with T-shaped ribs proposed in this study, the performance parameters are compared with the experimental data for the SAH with triangular ribs measured by Bekele et al. [8], as shown in Figures 19 and 20. The comparison was performed for the same cross-sectional area of the ribs and the same spacing between the ribs. The T-shaped rib used for the comparison is the reference rib presented in Table 2. Figure 19 shows the

variations of the Nusselt number (Nuo) and friction factor (*f* <sup>o</sup>) with Reynolds number. In the case of the friction factor, the T-shaped ribs show a slightly larger value than that with the triangular ribs in the entire Reynolds number range. However, the Nusselt number shows a large uniform improvement at all Reynolds numbers. Compared to the triangular ribs, the T-shaped ribs show Nusselt numbers that are 65% higher at Re = 4510 and 21.5% higher at Re = 22,600. On the other hand, a 35% higher pressure drop occurs at Re = 4510 and a 31% higher pressure drop occurs at Re = 22,600 than for the triangular ribs. Figure 20 shows the comparison of the performance factor PF between the two different ribs. The PF of the T-shaped ribs is 49.7% higher at Re = 4510, and 11% higher at Re = 22,600 than for the triangular ribs.

**Figure 17.** Variation of PF with h/e (rib area = 375 mm<sup>2</sup> ).

**Figure 18.** Nusselt number distributions on the heated surface between fifth and seventh ribs for different h/e. (**a**) h/e = 0.38; (**b**) h/e = 0.83; (**c**) h/e = 1.75.

Table 3 shows a comparison of PF ranges among various obstacle shapes developed to date. These obstacles were tested by previous experimental or numerical studies using the parameters of longitudinal pitch (P<sup>l</sup> /e), relative obstacle height (e/H), relative roughness height (e/D), relative roughness pitch (P/e), relative longitudinal pitch between the rows of winglets (P/H), and number of waves on the delta winglet (φ). The T-shaped ribs proposed in this study cause two recirculation regions formed in the height direction, as shown in Figures 8 and 13, due to their unique geometric shape. These multiple recirculation zones are expected to promote the production of turbulent kinetic energy and thus further enhance the turbulent heat transfer compared to the other obstacles in Table 3. **Figure 19.** Comparison of performance parameters between T-shaped and triangular ribs.

**Figure 20.** Comparison of PF between T-shaped and triangular ribs.

**Table 3.** Comparison of performance factor range obtained in the present study with those obtained in the previous experimental or numerical studies on obstacle shape.


### **6. Conclusions**

In this study, T-shaped ribs with a staggered arrangement are proposed to enhance the heat transfer in an SAH, and the effects of the geometric parameters of the ribs on the thermal and aerodynamic performance of the SAH were analyzed using RANS equations. The numerical results for the heat transfer rate and friction factor were validated by comparing with experimental data for the smooth duct and the SAH with triangular obstacles. In the case of the SAH with the obstacles, the average relative error of the average Nusselt number between the numerical and experimental results was less than 4.43% within the tested Reynolds number range, and that of the friction factor was less than 1.94%. In a parametric study using the ratios of the rib height to the channel height (h/H) and the rib width to the channel width (e/W), the performance factor (PF) was far more sensitive to h/H than to e/W, and the maximum PF was found at h/H = 0.5 for a fixed e/W. When h/e was varied by keeping the rib area constant, the variation of the friction factor with h/e was very small due to the fixed rib area, and thus the variation of the average Nusselt number was qualitatively similar to that of PF and their maximum values were commonly found at h/e = 0.83. The proposed T-shaped ribs showed superior performance to several other heat-transfer enhancement obstacles developed to date. For example, when compared with the triangular ribs, the T-shaped ribs showed a 49.7% higher PF at a Reynolds number of 4510, and an 11% higher PF at a Reynolds number of 22,600. To generalize the present results for staggered T-shaped obstacles, further studies are required for the SAH using different arrangements of T-shaped obstacles.

**Author Contributions:** S.-Y.A. presented the main idea of the T-shaped ribs; S.-Y.A. and K.-Y.K. contributed to the overall composition and writing of the manuscript; S.-Y.A. performed numerical analysis and analyzed the data; K.-Y.K. revised and finalized the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (No. 2019R1A2C1007657).

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

### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 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/).

## *Article* **Study on the Transient Characteristics of the Centrifugal Pump during the Startup Period with Assisted Valve**

#### **Qiao Li <sup>1</sup> , Xiang Ma <sup>2</sup> , Peng Wu 1,\* , Shuai Yang <sup>1</sup> , Bin Huang 3,\* and Dazhuan Wu 1,4**


Received: 19 August 2020; Accepted: 28 September 2020; Published: 2 October 2020

**Abstract:** The startup period, one of several transient operations in a centrifugal pump, takes note of some problems with these devices. Sometimes a transient high pressure and high flow rate over a very short period of time are required at the startup process. The pump' s dynamic response is delayed because of the rotational inertia of the pump and motor. Our research focuses on how to get a large flow in a short time when the pump cannot meet the requirements alone without a large power driver. To achieve a strong response in the startup process, a ball valve is installed downstream of the pump. The pump' s transient behavior during such transient operations is important and requires investigation. In this study, the external transient hydrodynamic performance and the internal flow of the pump during the transient startup period are studied by experiments and simulations. In order to find an appropriate matching method, different experiments were designed. The content and results of this paper are meaningful for performance prediction during the transient pump-valve startup period.

**Keywords:** transient characteristics; centrifugal pump; startup period; numerical simulation

#### **1. Introduction**

The transient operations of a pump such as the startup and shutdown periods occur in several applications of hydraulic systems. In a Navy application, the pump is used to launch weapons from submarines [1]; the pump starts up in a very short time to generate instantaneous pressure. In large pump stations, a sudden startup of the pumps will cause strong pressure fluctuations and power shocks.

The transient characteristics of pumps have been studied in recent years. Tsukamoto and Ohashi [2,3] studied the transient characteristics of a centrifugal pump, and the acceleration of the pump at the startup stage was very fast. They found that the hysteresis and transient pressure fluctuation around the vanes results in a difference between the dynamic and quasi-steady characteristics.

Lefebvre and Barker [1] conducted an experimental study on the acceleration and deceleration cycle of a mixed-flow pump. They found that the transient effect was notably on the pump performance. The quasi-steady assumptions are wrong for the performance prediction of impellers under highly transient conditions, such as at fast startup and stop periods. Saito [4] conducted an experimental study on the transient period of a pump. The results of the tests indicate that the locus of operating points on the head-capacity plane during the pump startup deviates from the system resistance curve to a great

degree as the mass of water in the pipeline becomes large. Dazin [5] and Wang [6–8] conducted startup and shutdown experiments for a pump. The transient characteristics were mainly influenced by the variation in transient velocity.

Wu [9,10] conducted simulations of the blades during startup periods. In two-dimensional simulations, he found that the sliding mesh method has a higher efficiency and stability than the dynamic mesh and dynamic reference frame methods. The startup period also influences the transient performance. Li [11,12] studied the transient characteristics of a centrifugal pump through experiments and numerical simulations during the startup period. He used the DSR (Dynamic Slip Region) method to determine the dynamic motion of the blades in the numerical simulation. The practicability of the DSR method in the dynamic impeller simulation was proved by the consistency between the experimental results and the numerical results. The instantaneous power shock is related to the joint action of a low flow rate and high acceleration in the initial stage of the startup, but the peak value is small. Due to the fluctuation in data, the observation is not very obvious.

In some applications, the pump needs to start up faster but the drive cannot provide sufficient starting acceleration, while in some applications the water hammer needs to be suppressed. A throttle valve can be installed on the pipeline, and the above-mentioned application needs can be met by adjusting the opening of the valve during startup.

The transient characteristics of valves have been studied extensively. Ming-Jyh Chern [13] used the particle tracking flow visualization (PTFV) method to conduct experimental research on the internal flow of a ball valve. He found that the valve' s flow coefficient is connected with the opening degree and has nothing to do with the velocity. Yang et al. [14] showed that the main pressure drop occurred along the valve throat because the circulation area diminished when the fluid flows through the valve throat. Wang [15] used the CFX software to conduct a three-dimensional dynamic analysis during the opening period of a butterfly valve. The transformation law of the transient flow coefficient was studied during the startup period. Cho [16] carried out an experiment and simulation to study a globe valve's force balance. The pressure distribution and stress distribution were studied.

Srikanth [17] studied industrial pneumatic valves by adopting dynamic meshing technology. Experiments and numerical simulations were carried out to research the startup period of a pneumatic valve. In the numerical simulation, the unstructured tetrahedral mesh was used to discrete the flow domain, the dynamic movement of valves was defined by the dynamic mesh method, and the standard k-epsilon model was implemented. He analyzed the effects of the opening degree and speed on the flow characteristics. The simulation and experimental results showed good similarity, thereby demonstrating the effectiveness. Quang Khai Nguyen [18] has conducted experiments to research the globe valve' s flow coefficient. He found that the R<sup>e</sup> and opening degree influences the flow coefficient.

The research objective of this article was to investigate a pump with an assistant valve. At present, there are few studies on the simultaneous startup of pump and valve. The key point of this paper is to study the simultaneous startup stage of a pump and assistant valve. In order to study the simultaneous startup stage of a pump and assistant valve, we designed many experiments and numerical simulations. Then, we analyzed the transient dynamic characteristics of the pump-valve system; the results of numerical simulations show the influence of the opening valve on the internal transient flow structure during the pump startup. Furthermore, we compare various startup modes and select the most suitable mode. In this paper, we focus on how to achieve a higher flow acceleration through a valve if the pump starts not fast enough.

#### **2. Experiment**

#### *2.1. Test Equipment and Method*

A schematic view of the experiment system is shown in Figures 1 and 2. The test system consists of a data acquisition system, a pump and control system, a valve, pipes, and data measurements. A ball valve is applied to change the dynamic characteristics by changing the valve close time. Two pipes are connected to the tank, and the pipe which connects the tank and water-tap is used to hold the water level in the tank. In order to guarantee that cavitation does not happen, the level is kept about 1.2 m above the pump inlet during this series of experiments.

**Figure 1.** Illustration of the experiment pump (vertical pipe pump). )

)

The data are acquired with an NI-4497. The data acquisition frequency is 5000 Hz. The pump comprises the motor, pump body, impeller, and diffuser. The parameters for the pump are listed in Table 1. The pump is controlled by a frequency conversion cabinet. The data measurements mostly contain pressure sensor, electromagnetic flowmeter, and photoelectric speed sensor. In the experiment, the pump starts up first; meanwhile, the valve remains closed for a period of time, *Tc*, and then the valve is rapidly full opened.

The installation locations of three pressure sensors are shown in Figure 2. The transient pressure at the inlet, outlet, and valve-behind section were measured.

The pump transported water from a tank to the atmosphere during the experiments. When the ball valve was fully opened, the flow rate was 0.0033 m<sup>3</sup> /s. Then, we changed the valve opening degree to adjust the flow rate to 0.0033, 0.00326, 0.00309, 0.0025, 0.00136, and 0 m<sup>3</sup> /s, respectively. The aim of

**Figure 2.** Diagram of the experiment (units are mm)**.**

the experiments was analyzing the steady performances at different flow rates. It took 9.1 s for the pump to reach 2900 rpm in this series of experiments.


**Table 1.** Specifications for the experimental pump.

The transient performances of the pump-valve system were studied by many experiments during the startup process. We started the pump first and then kept the valve closed until the pump' s rotational speed reached a certain value. The closed time for the valve was called Tc, and Tc was set as 0, 2.0, 3.2, 4.4, 5.9, 7.2, 8.0, 10.28, and 11.56 s, respectively. Then, at time point Tc, the valve was quickly opened in time To, and the data were collected until the pump reached the rated speed.

#### *2.2. Processing the Experiment Data*

A photoelectric speed sensor was used to measure the pump' s transient rotation speed. The photoelectric speed sensor outputs the pulse signal. The rotational speed used the pulses gathered in the time interval. The formula is expressed as:

$$m = \frac{60\text{N}}{\Delta t} \tag{1}$$

where *n* is the transient speed and *N* is the number of pulses in ∆*t*.

The transient head could be calculated through measuring the transient pressure of the inlet section and outlet section, and the transient head also contained a component which was used to accelerate the fluid in the startup period.

$$H = \frac{P\_o - P\_i}{\rho g} + \frac{Q^2}{2g} \left( \left(\frac{4}{\pi d\_o^2}\right)^2 - \left(\frac{4}{\pi d\_i^2}\right)^2 \right) + \left\{ L\_{eq} / (gA\_0) \right\} \frac{dQ}{dt} \tag{2}$$

where *Leq* is the equivalent length, *P*<sup>0</sup> is the outlet pressure, *P<sup>i</sup>* is the inlet pressure, and *A* is the cross-sectional area of the pipeline.

The dimensionless head and flow are defined as equaling:

$$\mathcal{C}\_H = \frac{P}{\rho u\_2^2 / 2} \tag{3}$$

$$\mathcal{C}\_{q} = \frac{\mathcal{Q}}{\pi d\_{2} b\_{2} u\_{2}} \tag{4}$$

where *u*<sup>2</sup> = π*d*2*N<sup>r</sup>* /60, *N<sup>r</sup>* is the current rotational speed of the pump, *b*<sup>2</sup> is the outlet width of the impeller, and d<sup>2</sup> is the diameter of the impeller.

Transient flow is an important parameter in the startup period. The flow is mostly measured by an electromagnetic flowmeter and an orifice flowmeter. When measuring the transient flow, the electromagnetic flowmeter has a greater delay. Figure 3 is an illustration of the orifice flowmeter. Bernoulli equation and the continuity equation between the high pressure and orifice are used to obtain the flow. Formula (5) shows the flow transition equation:

$$q\_V = \mathbb{C} \cdot Q\_V = \frac{\mathbb{C}}{\sqrt{1 - \beta^4}} \cdot \frac{\pi d^2}{4} \cdot \sqrt{\frac{2\Delta P}{\rho}} \tag{5}$$

where *C* is the outflow coefficient and β equals *d*/*D*. *β*

**Figure 3.** Illustration of the orifice flowmeter.

The outflow coefficient C is used to correct the pressure difference between the orifice and the low pressure. Yet, the orifice flowmeter does not consider the pressure loss of the orifice. In this paper, the local resistance loss is used to get the flow. Formulas (6) and (7) show the energy conservation equation and the local resistance equation. The local resistance coefficient is influenced by the orifice structure and the Reynolds number. Thus, it can be found that the pressure difference is a quadratic function of the flow:

$$\frac{P\_1}{\rho g} + \frac{u\_1^2}{2} = \frac{p\_2}{\rho g} + \frac{u\_2^2}{2} + h\_{f\bar{l}}\tag{6}$$

$$h\_{fl} = \xi \frac{l\_e}{d} \frac{u^2}{2} \tag{7}$$

2 where ξ is the resistance coefficient and *l<sup>e</sup>* is the equivalent length.

ξ Figure 4a is a schematic diagram. The distance between the low pressure and the orifice plate is 6D, and between the high pressure and the orifice it is 3D. The pressure measuring point for the local resistance flowmeter should be in a stable flow area, while the low pressure point in the orifice flowmeter should be in an unstable area. An electromagnetic flowmeter was installed downstream to get a precise steady flow. The orifice plates come in three sizes, and ultimately the first one was chosen because it causes the biggest pressure loss; therefore, in the small flow we can get a more definitive pressure difference. The pump is controlled by the inverter, then a series of flowrates can be found by changing the frequency. Through a series of steady state experiments, the relationship of pressure difference and flowrate can be found. Figure 5 shows the fitting curve and fitting formula.

(**b**)

**Figure 4.** Illustration of the local resistance flowmeter; (**a**) the schematic diagram, (**b**) three sizes for the orifice plate.

**Figure 5.** The fitting curve of the pressure difference and flowrate.

#### *2.3. Results and Discussion for the Experiment*

The rotational speed is shown in Figure 6. The startup time for the pump was 9.1 s. Figure 7 shows the changes in pressure during the startup period.

**Figure 6.** The rotational speed of pump measured by the photoelectric speed sensor.

**Figure 7.** *Cont*.

**Figure 7.** The inlet pressure, outlet pressure, and pressure behind the valve under different *T<sup>c</sup>* during the startup process; (**a**) *T<sup>c</sup>* = 0 s, (**b**) *T<sup>c</sup>* = 2.7 s, (**c**) *T<sup>c</sup>* = 4.2 s, (**d**); *T<sup>c</sup>* = 5.8 s, (**e**) *T<sup>c</sup>* = 7.8 s, (**f**) *T<sup>c</sup>* = 10.3 s.

The experimental process can be divided into five parts. Zone 1 means that the pump was stopped, zone 2 means that the pump was opened while the valve was closed, zone 3 means that the valve was opening, zone 4 means that the valve was fully opened, and zone 5 means that the pump system became stable. *T<sup>c</sup>* is the closing time for the valve. *T<sup>c</sup>* = 0 means that the valve was completely open during the pump startup period. When the valve was closed, the inlet pressure and pressure behind the valve remained unchanged, while the outlet pressure increased more rapidly than the operating condition at *T<sup>c</sup>* = 0. The outlet pressure continued to increase more rapidly. When the valve was in the process of opening, the outlet pressure had a sudden fall in 0.4 s. Meanwhile, the pressure behind the valve increased rapidly and a pressure fluctuation occurred at the inlet of the pump. When the valve was completely open, the pressure behind the valve kept growing slower. Finally, the pressures in the different working conditions reached the same value. The outlet pressure was greater than the pressure behind the valve. This deviation could be caused by the resistance between the valve and pipe. There was a jump in the rotational speed curve at the initial stage of the start that does not match with the pressure curve and head curve. In the initial stages of startup, the energy provided by the centrifugal pump was smaller than the resistance of pipe system, and therefore the response could not have been sensitive to the pressure.

In order to compare different working conditions, the dimensionless flow and head were analyzed. Figure 8 shows the changes in dimensionless flow during the startup period. *T<sup>o</sup>* is the startup time of the pump. The red curve is the startup period of the pump with a completely open valve. When the valve was closed, the dimensionless flow remained around 0. When the valve was opening, *C<sup>q</sup>* rapidly increased and the curve rapidly tended to the red completely open curve. When *T<sup>c</sup>* was less than *To*, the *C<sup>q</sup>* first increased rapidly and then had a slow increment with the same tendency compared with the completely open startup period. When *T<sup>c</sup>* was bigger than *To*, the *C<sup>q</sup>* quickly increased to the final value.

Figure 9 shows the changes in the non-dimensional head during the startup period. The red curve is the startup period of the pump with a completely open valve. In the pump startup period, when the valve remained closed *C<sup>h</sup>* was bigger than the completely open case. In the opening period of the valve, the *C<sup>h</sup>* had a sudden drop and then tended to the completely open case rapidly. Finally, *C<sup>h</sup>* in all those cases reached a certain value. When *T<sup>c</sup>* divided by *T<sup>o</sup>* was less than or equal to 0.46, the *C<sup>h</sup>* reached the maximum when the pump was fully started. When *T<sup>c</sup>* divided by *T<sup>o</sup>* was greater than 0.46, the *C<sup>h</sup>* reached the maximum at the point where the valve was going to open. When *T<sup>c</sup>* was less than *To*, the maximum of *C<sup>h</sup>* increased as *T<sup>c</sup>* increased. When *T<sup>c</sup>* was greater than *To*, the maximum of *C<sup>h</sup>* maintained a certain value. If we want to achieve a high initial *Ch*, then *T<sup>c</sup>* should be bigger than *To*.

**Figure 8.** Dimensionless flow for different startup periods.

**Figure 9.** Dimensionless head for different startup periods.

Considering the differential of the flow, Figure 10 shows the changes for different startup periods. The red curve is the startup period of the pump with a completely open valve. Figure 10b is the enlarged image of differential *Q* when *Tc*/*T<sup>o</sup>* = 0; when the valve was completely open during the startup period, the *dQ*/*dT* increased to the maximum of about 2.25 and then fell to 1.25 or so, and finally the *dQ*/*dT* fell to 0. From Figure 6, it can be found that the rotational speed has a higher acceleration at first and then transitions to a certain smaller acceleration. The change in speed reflects the change in *dQ*/*dT*. Figure 10a shows the *dQ*/*dT* for different startup periods. During the valve's opening period, the *dQ*/*dT*

quickly reaches the maximum and then falls to the red curve, and finally tends to 0. When *T<sup>c</sup>* is less than *To*, the maximum of *dQ*/*dT* increases when *T<sup>c</sup>* increases. When *T<sup>c</sup>* is bigger than *To*, the maximum of *dQ*/*dT* does not increase and even decreases a little when *T<sup>c</sup>* increases. Thus, in order to obtain the biggest flow during a given short period of time, the *T<sup>c</sup>* should be in the vicinity of 1.13 *To*, and there is no need to increase *T<sup>c</sup>* when *T<sup>c</sup>* is bigger than *To*.

**Figure 10.** Differential flow for different startup periods; (**a**) flow change rate (*dQ*/*dT*) for different startup periods; (**b**) flow change rate (*dQ*/*dT*) curve when *T<sup>c</sup>* is in the vicinity of 1.13 *To*.

− According to the above results, in order to find the most suitable *Tc*, a series of experiments when *T<sup>c</sup>* is in the vicinity of 1.13 *T<sup>o</sup>* have been performed. The results are shown in Figure 10b; it can be found that when *T<sup>c</sup>* = *To*−*Tv*, the biggest flow can be found during a given short period of time.

#### **3. Numerical Simulation**

#### *3.1. Computational Domain and Grid*

The overall view of the 3D model is shown in Figures 2 and 3. The numerical 3D model' s parameters were the same size as the real test system. The schematic view of numerical pump model and simplified valve model is shown in Figure 11. Table 2 shows the general view of the mesh of the pump-system model. The computational domain is discretized by an unstructured tetrahedral cell. The grid independence test is carried out, and Figure 12 shows the results. In the calculation, the scaled residuals of continuity equation and momentum equations are reduced to a magnitude below 10−<sup>5</sup> , and the scaled residuals of the turbulence kinetic energy (*k*) and dissipation rate (ξ) are reduced to a magnitude below 10−<sup>4</sup> . The y+ value near the critical surfaces is from 0 to 30 and the average value is within 50. The difference in head between the 1.98 <sup>×</sup> <sup>10</sup><sup>6</sup> grid number and the 2.71 <sup>×</sup> <sup>10</sup><sup>6</sup> grid number is 0.011%, and the difference in flow is 0.026%. As a result, the total number of grids is ensured to be 1.98 <sup>×</sup> <sup>10</sup><sup>6</sup> . <sup>−</sup> *ξ* − <sup>−</sup> *ξ* −

**Figure 11.** The numerical model for the pump and ball valve.

**Table 2.** Overview of the mesh in the numerical model.

Steady simulations were performed firstly. The *Q-H* curve between the experiment and steady simulation is shown in Figure 13. The strong similarity in the results verifies the accuracy of the 3D model and grid quality. The flow in the follow-up experiment and transient simulation was around 0.003 m<sup>3</sup> /s.

**Figure 13.** Comparison of the experimental and steady numerical simulation.

#### *3.2. Boundary Conditions of the Transient Simulation*

2900

 

The startup period of *T<sup>c</sup>* = 8 s was simulated. Li [11] demonstrated the applicability of the detached eddy simulation (DES) model and the dynamic slip region (DSR) method in a transient simulation of the pump's startup periods. For the DES (detached eddy simulation) model, in the attached boundary layer, the RANS equation is applied; the large eddy simulation (LES) model is applied in other regions. The LES is more precise and the RANS is more efficient, and the DES model integrates the advantages. Using the DES model can save computational resources and computing time to some extent.

The DSR (dynamic slip region) method segments the flow domain into dynamic and static regions. Then, the dynamic region and stationary region are discretized, respectively, and then a pair of interfaces is used to transfer and exchange data to the two individual parts. The interfaces could be set up in the fluent.

A mathematical curve is used to define the speed. It was achieved with measurements from the experiment. Figure 6 shows the rotational speed curve, and formula 8 (the unit is r/min) shows the definition of speed.

$$\text{speed} = \begin{cases} 2747.65468 \times t & t \le 0.1 \\ 317.495 \times t^2 - 41.039 \times t + 275.315 & 0.1 < t \le 0.85 \\ 294.566 \times t + 219.44 & 0.85 < t \le 9.1 \\ 2900 & t > 9.1 \end{cases} \tag{8}$$

9.1

The valve's opening time is 0.4 s, and valve's rotational speed is assumed to be constant. Thus, the valve's rotational speed can be defined as follows (the unit is degree/s): 0 225 8 8 8.4

$$
\omega = \begin{cases} 0 & t < 8 \\ 225 & 8 \le t \le 8.4 \\ 0 & t > 8.4 \end{cases} \tag{9}
$$

In order to ensure both accuracy and efficiency, we set the time steps to vary with pump's rotational speed. Thus, time step is also a value that varies with time. Formula 10 (the unit is s) shows the definition of the time step. At each time step, the blade does not rotate more than one degree when t ≤ 0.85 s, the blade does not rotate more than four degrees when 0.85 s < t ≤ 9.1 s, and then the blade does not rotate more than six degrees. ≤ ≤ 0.000607 0.85

$$
\Delta t = \begin{cases}
0.000607 & t < 0.85 \\
\frac{\pi}{1388.1105 \times t + 1034.091} & 0.85 < t \le 9.1 \\
0.000345 & t > 9.1
\end{cases} \tag{10}
$$

The boundary conditions of the pressure inlet and pressure outlet were chosen. Second-order implicit was applied for the time-dependent term format. The SIMPLEC algorithm was used to compute the pressure–velocity coupling. The dispersion of convective terms was calculated by the second-order upwind scheme, and the numerical under-relaxation and diffusion terms were calculated by central difference schemes.

There is an assumption that the water was incompressible in the numerical simulations at the startup period and there was no cavitation. The numerical simulation utilized Ansys Fluent.

#### *3.3. Comparison of the Simulation and Experiment*

Figures 14–16 show the comparison of the transient simulation and the experiment. The experiment was recorded at *T<sup>c</sup>* = 7.8 s, while the simulation was made at *T<sup>c</sup>* = 8 s. The simulated pressure is probably about 2000 Pa higher than the experiment at the initial stage; this is caused by the measurement deviation in the water level in the tank. From the head and flow pictures, we can see that before the opening of the valve the simulation and experiment show a strong agreement, and the difference between the simulation and the experiment are mainly caused by a deviation of 0.2 s. In this experiment, the ball valve is controlled manually and the deviation is inevitable. Yet, the similarity between the numerical simulation and the experiment can prove the accuracy of the simulationn.

**Figure 14.** Comparison of the experiment and simulation pressures.

**Figure 15.** Comparison of the experiment and simulation heads.

**Figure 16.** Comparison of the experiment and simulation flows.

#### *3.4. Results and Discussion of the Transient Simulation*

The evolution of the velocity field at the opening period of valve is shown in Figure 17. A high-velocity zone appears at nearby of interface between valve and pipe when valve is opening. With the increase of opening degree, this high-velocity region decreases. The high-velocity region results in turbulence in the vicinity region, but there is little effect on the internal flow of pump because the valve's inlet flow is generally uniform. The change in resistance has a greater impact on the transient internal flow of pump. The sudden drop in head and pressure occurred because the pump works on the domain behind the valve.

**Figure 17.** The velocity evolution of the valve during the valve opening process: (**a**) *t* = 8 s, (**b**) *t* = 8.1 s, (**c**) *t* = 8.2 s, (**d**) *t* = 8.3 s, (**e**) *t* = 8.4 s.

The internal velocity evolution of pump during the opening period of valve is shown in Figure 18. In Figure 18a, it can be found that there is a clear vortex in volute tongue and trailing edge of blades. In Figure 18b,c, it can be found that the turbulence flow occurs in the leading edge of the blades. In Figure 18e, the pump's internal flow structure changes back to a representative transient flow structure when the valve is fully open. During the valve's opening period, the vortexes in the trailing edge of the blades and the volute tongue gradually disappears. In the valve's opening period, the performance of the pump has been improved. When the valve is closed, the pump's performance is weakened. This might be the reason why the maximum of flow change rate (*dQ*/*dT)* does not increase and even decreases a little, while *T<sup>c</sup>* increases when *T<sup>c</sup>* is bigger than *To*.

**Figure 18.** The internal velocity evolution of the pump during the valve opening process: (**a**) *t* = 8 s, (**b**) *t* = 8.1 s, (**c**) *t* = 8.2 s, (**d**) *t* = 8.3 s, (**e**) *t* = 8.4 s.

#### **4. Conclusions**

The research emphasis was a pump-valve system (centrifugal pump and ball valve). Many experiments and numerical simulations were carried out to study the transient performances of the pump-valve system. The study of transient performance mainly focused on how to get a large fluid acceleration.

In the experiments, the local resistance loss was used to measure the transient flow more accurately. The results of experiments show that a high and rapid transient response can be achieved by controlling the assistant valve during the pump's startup period. When *T<sup>c</sup>* > *To*, the maximum of *dQ* and *C<sup>h</sup>*

will not increase. Thus, the most suitable startup method for the pump and valve will be around *T<sup>c</sup>* = *To*. Based on the experimental results, when *T<sup>c</sup>* = *T<sup>o</sup>* − *T<sup>v</sup>* the biggest flow can be obtained during a given short period of time. Thus, the suggested suitable startup method for the pump and valve is *T<sup>c</sup>* = *T<sup>o</sup>* − *Tv*, which can get the biggest flow and a lower max pressure.

The most suitable *T<sup>c</sup>* is noted by considering the transient response's sensitivity and stability synthetically. When the most suitable *T<sup>c</sup>* is determined by the experiments, a numerical simulation is carried out to study the internal flow. The dynamic movement of the impeller and the valve is defined by the DSR method. The transient numerical simulation uses the DES model. The simulated results and experimental results show a good agreement. The simulation predicts the pressure and head volatility because of the valve's sudden opening. During the valve's opening period, the impeller-volute interaction plays a decreasing role in the pump's performance, while the vortex's revolution plays an increasing role. During the opening period of valve, the change in the valve resistance influences the internal flow of the pump.

This paper only investigated the consequence of *Tc*. Others are involved in the matching method, such as the valve degree. The key point of future work will be to study these factors synthetically and find the best matching mode.

**Author Contributions:** Conceptualization, Q.L. and P.W.; methodology, Q.L.; software, Q.L. and S.Y.; validation, Q.L., D.W., and X.M.; formal analysis, Q.L.; investigation, Q.L.; resources, D.W. and X.M.; data curation, Q.L.; writing—original draft preparation, Q.L.; writing—review and editing, P.W. and B.H.; supervision, P.W. and B.H.; project administration, P.W. and X.M.; funding acquisition, D.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Natural Science Foundation of China (No. 51839010, No. 52076186).

**Acknowledgments:** The author sincerely thanks the support of the State Key Laboratory of Fluid Power Transmission and Control of Zhejiang University.

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

#### **Nomenclature**


#### **References**


© 2020 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/).

### *Article* **Numerical Simulation of Axial Vortex in a Centrifugal Pump as Turbine with S-Blade Impeller**

### **Xiaohui Wang \*, Kailin Kuang, Zanxiu Wu and Junhu Yang**

School of Energy and Power Engineering, Lanzhou University of Technology, Lanzhou 730050, China; kuangkailin@163.com (K.K.); wzx19960512@163.com (Z.W.); lzyangjh@lut.cn (J.Y.)

**\*** Correspondence: wangxh@lut.edu.cn

Received: 18 August 2020; Accepted: 11 September 2020; Published: 20 September 2020

**Abstract:** Pump as turbines (PATs) are widely applied for recovering the dissipated energy of high-pressure fluids in several hydraulic energy resources. When a centrifugal pump operates as turbine, the large axial vortex occurs usually within the impeller flow passages. In view of the structure and evolution of the vortex, and its effect on pressure fluctuation and energy conversion of the machine, a PAT with specific-speed 9.1 was analyzed based on detached eddy simulation (DES), and the results showed that vortices generated at the impeller inlet region, and the size and position of detected vortices, were fixed as the impeller rotated. However, the swirling strength of vortex cores changed periodically with double rotational frequency. The influence of vortices on pressure fluctuation of PAT was relatively obvious, deteriorating the operating stability of the machine evidently. In addition, the power loss near impeller inlet region was obviously heavy as the impact of large axial vortices, which was much more serious in low flow rate conditions. The results are helpful to realize the flow field of PAT and are instructive for blade optimization design.

**Keywords:** energy recovery; pump as turbine; vortex; hydraulic losses; pressure fluctuation

#### **1. Introduction**

In recent decades, pump as turbine (PAT) has drawn increasing attention in energy recovery systems where a high-pressure water source exists. With a pump operating as turbine, the direction of flow and rotation are opposite. For pump, energy is supplied to the fluid via a rotating shaft, as shown in Figure 1a, it is a energy absorbing device. For PAT, energy is extracted from the fluid and output via the rotating shaft, as shown in Figure 1b, it is a energy producing device. Compared with a conventional hydraulic turbine, PAT is simple, inexpensive, easy to maintain, readily available worldwide, and has a short capital payback period. It is an attractive solution for micro-hydro power with capacity below 100 kW [1]. It would be economical to use PAT, recovering the dissipated energy of high-pressure fluids in several hydraulic energy resources, such as water distribution network (WDN) [2–5], sea water reverse osmosis (SWRO) [6], chemical processes [7], nature falls [8], etc.

Although there is a wide application of PAT, the selection of a proper pump operating as turbine is particularly challenging. Many selection techniques have been published so far, while researchers have tested their models on few pumps and recorded deviations in the order of ±10~20% [9]. PAT may not have optimum or favorable flow behavior since pumps are usually not designed for turbine operation. The mismatch between turbine flow parameters and pump geometry may affect the stability of flow performance [10]. In addition, PATs have poor part-load performances [11,12]. Many researchers have presented the optimization of the turbine mode performance for overcoming these challenges. All mentioned issues above require a detailed understanding of the internal flow mechanism of PAT, which is important and imperative to predict the performances of PAT, as well as improve its efficiency and operating stability.

**Figure 1.** Directions of flow and rotation in pump and pump as turbine (PAT). (**a**) Pump; (**b**) PAT.

Computational fluid dynamics (CFDs) were adopted extensively in many earlier research studies to investigate the flow behavior of PAT. A detailed analysis of the turbulence flow structure of a pump and its reverse mode were performed by Pascoa et al. [13]. Singh and Nestmann [14] revealed the wakes and the corresponding losses (flow separation) at the inlet and exit of a PAT impeller. Yang et al. [15] discussed the velocity distribution and hydraulic losses of PAT with different blade wrap angles. Ardizzon and Pavesi [11] researched the effect of relative through-flow and eddy vortex on flow behavior in PAT impellers and established the optimum incidence angle in outward- and inward-flow impellers. Zhang et al. [16] presented a numerical simulation study that the reverse flow causes a great deal of vortex in impellers. Stanbli et al. [17] studied the instability of PAT during start-up process using numerical-based method. Zobeiri et al. [18] investigated the rotor–stator interactions in turbine mode of a pump and presented the pressure fluctuation in stator flow channels. Singh and Nestmann [9] analyzed the flow condition in different flow zones of PAT and concluded the hydraulic loss of each flow zone. Simão et al. [19,20] investigated the hydrodynamic flow behavior of centrifugal PAT to better understand the energy recovery system behavior and to reach the best efficiency operation conditions. Additionally, collaborative design of rotor and stator of PAT have been concerned to improve its efficiency recently [21,22].

– The PATs have poor hydraulic performances usually as the pump manufactures do not pay attention to the performances of a pump in reverse operation. As a consequence, low efficiency and instability have been found generally due to the poor flow conditions of vortices, secondary flow, and pressure fluctuation. In this research, the flow behavior of PAT was simulated by the CFD method. The structure and evolution of the large axial vortex in impeller channels were revealed, and its impact on pressure fluctuation and power losses was observed. The results can be expected to be a support for the optimization design of PAT.

#### **2. Theoretical Model of Vortex**

For PAT, the slip phenomenon occurs inevitably in flow passages caused by finite blades. Due to the finite blades with certain thickness of PAT, the fluid in the flow passages is guided weakly, and subsequently a slip velocity ∆*c<sup>u</sup>* is generated, as shown in Figure 2 (where *u* is the peripheral velocity of impeller, *w* is relative velocity, *c* is absolute velocity, *c<sup>m</sup>* and *c<sup>u</sup>* are the meridian and peripheral components of absolute velocity, respectively, ∆*c<sup>u</sup>* is slip velocity, β is relative flow angle and β*<sup>b</sup>* is the blade angle, the subscript ∞ represents infinite blades). As a consequence, the large axial vortices induced reasonably. In part-load operation, this phenomenon is much more serious on account of the non-optimum incoming flow.

*Δc*

*β* is the blade angle, the subscript ∞ represents infinite blade

**Figure 2.** The model of the large axial vortex in PAT.

In recent decades, various vortex identification methods, including closed or spiral path lines, minimum local pressure, vorticity magnitude, etc., have been used to interpret vortical structures in instantaneous flow fields [23]. The vorticity magnitude is widely applied in qualifying the intensity of vortices, and the swirling strength has been adopted—usually to qualify the vorticity magnitude.

For a random element of vortex flow, the velocity gradient tensor *dij* can be described as

$$
\begin{bmatrix} d\_{ij} \end{bmatrix} = \begin{bmatrix} v\_r v\_{cr} v\_{ci} \end{bmatrix} \begin{bmatrix} \lambda\_r & 0 & 0 \\ 0 & \lambda\_{cr} & \lambda\_{ci} \\ 0 & -\lambda\_{ci} & \lambda\_{cr} \end{bmatrix} \begin{bmatrix} v\_r v\_{cr} v\_{ci} \end{bmatrix}^{-1} \tag{1}
$$

*λ*

—

*Δc β*

 where *v<sup>r</sup>* , *vcr* , and *vci* represent the axial, radial, and tangential components of the element velocity, respectively. The velocity gradient tensor *dij* exists one real eigenvalue λ*<sup>r</sup>* and two conjugated complex eigenvalues λ*cr* ± λ*ci* . The swirling strength is the imaginary part of the complex eigenvalues of the velocity gradient tensor, λ*ci* ; it is positive if and only if the discriminant is positive and its value represents the strength of swirling motion around local centers. The greater the absolute value of the swirling strength, the stronger the internal circulation of fluid.

#### *λ λ* **3. Numerical Simulation**

#### *λ 3.1. Numerical Method*

A modified PAT with specific speed (*nQ*1/<sup>2</sup> /*H*3/<sup>4</sup> , where *n* is rotational speed, *Q* is flow rate, and *H* is head) 9.1 was selected for numerical simulation, and the hydraulic parameters were 50 m for head and 50 m<sup>3</sup> /h for flow rate with the rotational speed 1500 r/min. The flow zones consisted of volute, impeller, and draft tube as shown in Figure 3. The main geometrical parameters were shown in Table 1.

**Figure 3.** Structure of modified PAT. (**a**) Structure of selected PAT; (**b**) modified impeller of PAT.



– The numerical simulation was performed by means of the Navier–Stokes equation with an appropriate turbulence model. The Reynolds Averaged Navier-Stokes (RANS) turbulence model is not appropriate for the unsteady flow prediction, while a fully-resolved Large Eddy Simulation (LES) is almost unfeasible nowadays [24]. Recently, the Detached Eddy Simulation (DES) showed the superiority of the prediction of unsteady flow phenomenon in studies by Magnoli and Schilling [25]. DES can be described as a hybrid RANS-LES turbulence modeling approach and can be applied in a numerical simulation of rotor–stator interaction, inter-blade vortices and vortex rope successfully. It is acting as a Sub-Grid-Scale (SGS) model of LES in regions where the grid resolution is fine enough to resolve turbulent structures, while in other regions the model is used as a pure RANS model [26]. It features the advantages of a less refined grid near the wall, as well as the memory requirements of a computer. DES can be explicitly presented in the Spalart–Allmaras (SA) *k*–ε model or Shear Stress Transport (SST) model.

– In the present work, the CFX (17.0, ANSYS, Pittsburgh, PA, USA, 2016) was adopted for the solution of 3D Navier–Stokes equations due to its characteristics of robust and fast convergence [15]. Steady simulations were achieved using the RNG *k*–ε model and the results were applied as the initial value of transient analysis, the SA-DES turbulence model was applied for transient simulations.

The one equation SA-DES model can be described as

$$\frac{\partial \overline{\boldsymbol{v}}}{\partial t} + \boldsymbol{u} \cdot \nabla \overline{\boldsymbol{v}} = \frac{1}{\sigma \text{Re}} [\boldsymbol{\nabla} \cdot ((\boldsymbol{v} + \overline{\boldsymbol{v}}) \nabla \overline{\boldsymbol{v}}) + c\_{b2} [\boldsymbol{\nabla} \overline{\boldsymbol{v}}]^2] + c\_{b1} \overline{\boldsymbol{S}} \overline{\boldsymbol{v}} - \frac{c\_{w1} f\_{\overline{w}}}{\text{Re}} (\frac{\overline{\boldsymbol{v}}}{d\_{\text{DES}}})^2 \tag{2}$$

–

−

–

where <sup>e</sup> *<sup>v</sup>* is a destruction term for the eddy viscosity, which is proportional to (<sup>e</sup> *v*/*d*) 2 , where *d* is the distance to the closest wall, *Re* is Reynolds number. The second and last terms on the right side of the equation are the product term and destruction term, respectively. When balanced with the production term, the eddy viscosity is adjusted to scale with the local deformation rate *S* and *d*: e *v* ∝ *Sd* 2 . In the Smagorinsky model, the sub-grid-scale (SGS) eddy viscosity scales with *S* and the grid spacing <sup>∆</sup>:µ*SGS* <sup>∝</sup> *<sup>S</sup>*<sup>∆</sup> 2 . Thus, the SA model turns into the SGS model when *d* is replaced by a length proportional to ∆. grid spacing Δ: Δ.

–*ε*

If we replace *d* in the SA destruction term with *d* e , described as

$$
\tilde{d} = \min(d, \mathbb{C}\_{DES} \Delta) \tag{3}
$$

then the model is an SA turbulence model when *<sup>d</sup>* <sup>≪</sup> <sup>∆</sup>, while an SGS model when *<sup>d</sup>* <sup>≫</sup> <sup>∆</sup>. 

Since the SA-DES model does not require any wall functions, the mesh that is close to the wall surface must be designed to accurately predict the hydrodynamic force; hence, the high-aspect-ratio cells near the wall have been generated [27]. The mesh of the fluid domain was generated using ICEM-CFD (17.0, ANSYS, Pittsburgh, PA, USA, 2016) as its advantage of a well-adapted and efficient hexahedral grid was applied for meshing, as shown in Figure 4. The length of the inlet and outlet pipes was extended to eliminate the influence of back flow. The grid convergence and grid independence tests were performed, and the results showed that the head–flow rate curve became stable as the elements of mesh over 3,612,548, as shown in Figure 5. Therefore, the final element numbers of the volute, impeller, and draft tube were 865,260, 3,007,154, and 243,200, respectively. –

**Figure 4.** Computational domain and mesh. (**a**) Computational domain; (**b**) mesh.

**Figure 5.** Grid independent test.

μ

−5

The boundary condition of the inlet was the total pressure with an initial value of 0.5 MPa, and the outlet was the flow rate with an initial value of 50 m<sup>3</sup> /h for design condition; the rotational speed of the impeller was fixed with 1500 r/min. The fluid was the normal water with a temperature of 20 ◦C, all the wall surfaces were adiabatic, and the roughness was set to 50 µm. To obtain reasonable results, the proper selection of time steps is of great importance. It is suggested that time steps for a runner rotation of 0.5–5 ◦ could provide useful information for the flow field under transients [28]. Hence, the time steps in this study were 3.3 × 10 −4 s, corresponding to 3 ◦ of the impeller rotational angle. The max coefficient loop of convergence control was 40, and the residual target of the convergence criteria was 10−<sup>5</sup> . The total time of the duration data was 0.4 s corresponding to 10 rotor revolutions and the last four revolution data were analyzed. μ – − −5

For revealing the vortex structure and pressure fluctuation in PAT flow channels, 17 monitoring points were set in the middle plane of PAT, as shown in Figure 6. Point 1 was set in the gap between volute and impeller, point 2 and 3 next to the inlet and outlet of impeller, respectively, point 4 was in center of impeller outlet, points 5, 7, 9 were on the suction profile of the short blade, while 6, 8, 10 were on the pressure side, points 11, 13, 15 were on the suction profile of the long blade, while 12, 14, 16 were on the pressure profile, point 17 was set in the flow channel.

**Figure 6.** Monitoring sites of PAT.

#### *3.2. Verification of Numerical Method*

In this section, the numerical method was validated by the experimental results. A test rig was established for the hydraulic performance experiment of PAT. The test rig was composed by water supply, PAT, and energy dissipation sections as shown in Figure 7. A feed pump was installed to provide the head and flow rate for PAT. A magnetic power brake was equipped to balance the output power, and a loop control system was used to adjust the torque of output shaft. A flow meter was equipped at the inlet pipe of PAT for measuring the flow rate, and two pressure transducers were installed at the PAT inlet and outlet for measuring pressure. For measuring the torque and rotational speed of PAT, a torque meter was set at shaft. The head, flow rate, power and efficiency of PAT could be obtained after all parameters were measured.

(**b**)

**Figure 7.** Experimental equipment of PAT. (**a**) Schematic diagram of experiment; (**b**) test rig.

The selected PAT was tested and the hydraulic performance curves by experimental and numerical methods were illustrated in Figure 8. It can be found that the numerical head is in good coincidence with the experimental results. In consequence, it is reasonable to believe that the employed numerical method is accurate, and it can be applied in performance predictions of PAT.

**Figure 8.** Comparison between experimental and numerical results.

#### **4. Results and Discussion**

#### *4.1. Vortex Information in Flow Channels*

For PAT, large axial vortices were derived in impeller flow channels even at the best efficiency point (BEP), as shown in Figure 9. It can be seen that the streamline was disordered near the suction surfaces of the impeller inlet where the large axial vortices were induced. It can be clearly observed

that the vortices were more legible on short blade surfaces. Furthermore, the size and position of detected vortices were invariant with rotor rotating, and apparently, these were stable.

**Figure 9.** Streamline and swirling strength contour of impeller.

−1 −1 −1 The swirling strength contour of detected vortices during one rotating cycle was displayed in Figure 9. It can be seen that although the vortex size and position were stable, the swirling strength changed with rotation of the impeller. In the first half of the rotating cycle, the swirling strength was minimum (0~50 s −1 ) at 0.122011 s, intensified to maximum (450~500 s −1 ) at 0.133891 s, and weakened to minimum (0~50 s −1 ) again at 0.142141 s. The revolution of the vortex swirling strength in the second half of the cycle was the same as the first half. The vortices information was extremely similar at 0.122011 s and 0.142141 s, 0.127951 s and 0.148081 s, as well as 0.133891 s and 0.154021 s. It can be found that the time steps of each working point are 0.2 s approximately for the three groups (group a, b, and c, as shown in Figure 9), that is half the time of per rotating cycle (0.4 s) for PAT. In other words, the swirling strength of the vortex develops periodically with two times the rotating frequency.

#### *4.2. Pressure Fluctuation of Vortex*

The pressure fluctuation could be produced due to the vortices with twofold rotating frequency. In order to reveal the pressure fluctuation characteristics of PAT, transient numerical simulation was performed, and results of the 17 monitoring points are given in Figure 10. Where the vertical coordinates *Cp* is pressure coefficient, it can be described as

$$\mathcal{C}\_p = \frac{p\_i - \overline{p}}{0.5\rho u\_1^2} \tag{4}$$

where *p<sup>i</sup>* denotes the transient pressure of monitoring point (Pa), *p* is the average pressure (Pa), ρ is the density of fluid (kg/m<sup>3</sup> ), and *u*<sup>1</sup> is the peripheral velocity of impeller inlet (m/s).

*ρ*

**Figure 10.** Pressure fluctuation coefficient with time. (**a**) Points 1, 2, 3; (**b**)Points 5, 7, 9; (**c**)Points 4, 6, 8; (**d**)Points 11, 13, 15; (**e**)Points 12, 14, 16; (**f**)Points 4, 17.

– – The pressure of point 1 fluctuated 10 times in a rotating cycle visibly, which was caused by the blade–volute interaction; as the point was set in the gap between the rotor and volute, it was not related to the axial vortex obviously. Point 4 could not be related to the axial vortex as well because it was set in the draft tube that was far away from the vortex regions. The pressure fluctuation at other points showed that the leading periodical impulse was related to the rotor–stator interaction. However, it was no reason to neglect the correlation between the subordinate periodical impulse and axial vortex.

Figure 11 showed the pressure fluctuation images with frequency range, which was received by fast Fourier transform (FFT) from Figure 9. For easily understanding, the horizontal axis was the ratio of frequency (*f*/*fn*), where *fn* denoted the rotating frequency of the rotor.

It can be seen that the leading pressure fluctuation of the monitoring points occurred at 1 *f*/*fn*, 5 *f*/*fn*, 10 *f*/*fn*, and 20 *f*/*fn*. Obviously, this related to the rotor–stator interaction. Pressure fluctuation at 1 *f*/*fn* caused by the rotor–tongue interaction, at 5 *f*/*fn*, 10 *f*/*fn*, 20 *f*/*fn* caused by blade–tongue interaction (the impeller equipped with 5 long blades and 5 short blades). Therefore, the main factor of pressure fluctuation was the rotor–stator interaction.

However, the subordinate pressure fluctuation at 2 *f*/*fn* was found in points 2, 5, 6, 11, 12, and 17. It was obvious, especially in points 5 and 11 as marked with dashed circle in Figure 10, as mentioned earlier, that the large axial vortices were derived in these regions usually. Consequently, it was reasonable to declare that the subordinate pressure fluctuation was related to axial vortices in impeller channels, which deteriorated the operating stability of the machine evidently.

–

*ρ*

–

**Figure 11.** Pressure fluctuation coefficient with frequency. (**a**) Points 1, 2, 3; (**b**)Points 5, 7, 9; (**c**)Points 4, 6, 8; (**d**)Points 11, 13, 15; (**e**)Points 12, 14, 16; (**f**)Points 4, 17.

#### *4.3. Power Losses Caused by Vortex*

– – – The large axial vortices provide subordinate contribution to the pressure fluctuation of PAT. More importantly, this might cause entropy generation in the flow field, and therefore, the power loss produced inevitably. In this section, the power losses caused by axial vortices were analyzed.

– Flow distortion have been detected in the impeller that was caused by axial vortices, where a wake region has been found near the impeller inlet, it was significant especially for 0.6 *Q<sup>d</sup>* and 1.0 *Q<sup>d</sup>* (*Q<sup>d</sup>* is design flow rate), as shown in Figure 12. As a consequence, a low-pressure zone appeared near the impeller inlet, and the relative velocity no longer distributed alongside the blade surfaces. In order to reveal the effect of axial vortices on performance characteristics of PAT, six monitoring cylindrical surfaces in the impeller were created as shown in Figure 13. Figure 14 was the distribution of the average relative velocity (radial component) in the impeller, and Figure 15 was the average pressure at each cylindrical surface. It can be seen that the average relative velocity (radial component) and pressure curves decreased gradually along the flow direction in the impeller channels for 1.6 *Qd*; however, a local decline of the curves appeared at surface 1 and 2 for 0.6 *Q<sup>d</sup>* and 1.0 *Qd*. As shown in Figure 12, the streamline of 1.6 *Q<sup>d</sup>* was uniform, and very tiny axial vortices were detected in the flow channels. However, large axial vortices can be found near the impeller inlet at 0.6 *Q<sup>d</sup>* and 1.0 *Qd*; this was much more serious for low flow rates. It can be seen from Figure 12 that the region from surface 1 to surface 3 was worse affected by axial vortices for the low flow rates that reasonably responded to the local decline of the average relative velocity and pressure.

**Figure 12.** Streamline at radial plane of impeller. (**a**) 0.6 *Q<sup>d</sup>* ; (**b**) 1.0 *Q<sup>d</sup>* ; (**c**) 1.6 *Q<sup>d</sup>* .

**Figure 13.** Monitoring surfaces of the impeller.

**Figure 14.** Average relative velocity of monitoring surfaces.

**Figure 15.** Average pressure of monitoring surfaces.

 (

(

( ( 

*ω*

(

*ρ*

As the axial vortices were generated, power losses were raised inevitably. To study the power losses caused by vortices, the flow domain in the impeller was divided into six zones (Figure 16), and power losses of each zone were calculated.

**Figure 16.** Zones of the impeller.

For any zone *i*, when the boundary condition with a pressure inlet and flow rate outlet are given, while rotating speed is fixed, the theoretical power (fluid power) and actual power (shaft power of PAT) can be obtained by numerical simulation. The theoretical power can be described as

$$p'\_{(i,i+1)} = \rho \text{gQH}\_{(i,i+1)}\tag{5}$$

*ρ* where ρ is the fluid density, *g* is the gravitational acceleration, *H*(*i"i*+1) is the fluid head of zone *i*, and *Q* is the flow rate. The shaft power of PAT is

$$
\sigma\_{(i,i+1)} = \mathcal{M}\_{(i,i+1)} \cdot \boldsymbol{\omega} \tag{6}
$$

*ω* where *M*(*i,i*+1) is the torque of zone *i*, while ω is the angular speed of the impeller. Then, the relative power losses of zone *i* are

( (

$$f\_{(i,i+1)} = 1 - \frac{p\_{(i,i+1)}}{p'\_{(i,i+1)}} \tag{7}$$

 ( ( As the numerical simulation did not consider the leakage and frictional losses of PAT, Equation (7) can be considered as relative power losses of zone *i.*

Figure 17 presented the power losses of each zone for 0.6 *Qd*, 1.0 *Qd*, 1.6 *Qd*, respectively. It can be seen that the power losses of zone 1 and 2 (the inlet region of the impeller) were higher distinctly than zone 3 and 4. As mentioned earlier, the large axial vortices occurred in this region usually, and caused the reduction in energy conversion of PAT. What calls for special attention was that the power losses of zone 5 and 6 were higher than zone 3 and 4 as well, which was related to the vortices in the draft tube to a great extent, and it deserved further research in the future.

It can be concluded also that the power losses were heaved significantly in low flow rates, which was related to the large axial vortices. As mentioned above, large axial vortices can be found near impeller inlets usually and are much more serious for low flow rates. Thus, it was believed that power losses would be induced by the large axial vortices within flow passages, and it should be considered in the design and optimization process, especially in low flow rates.

**Figure 17.** Hydraulic losses of different zones.

#### **5. Conclusions**

In this study, the flow behavior of a centrifugal PAT with specific speed 9.1 was researched by a verified CFD method. The large axial vortices were derived in impeller flow channels due to the slip and poor match between flow and blades. The size and position of the vortices were stable apparently. However, the swirling strength developed periodically with 2*fn* (*fn* is the rotating frequency of PAT).

The pressure fluctuation can be found in PAT. The leading pressure fluctuation was caused by the rotor–stator interaction, while the subordinate fluctuation was related to axial vortices in impeller channels, which deteriorated the operating stability of the machine evidently.

– The power losses were induced by the large axial vortices in the impeller flow channels, and this phenomenon was much more serious in the low flow rate operation. This should be considered in the design and optimization process, especially in low flow rates.

Nevertheless, the feature of large axial vortices within PAT impeller channels deserved further research in detail based on more PATs. The influence of the rotation, geometry of the blades, entropy variation of the large axial vortex and efficiency of the PAT should be discussed deeply, especially in the pump–turbine transition processes.

– **Author Contributions:** X.W.: Conceptualization, Methodology, Formal Analysis, Investigation, Writing. K.K.: Software, Data Curation. Z.W.: Resources, Data Curation. J.Y.: Validation. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Nature Science Foundation of China (51569013), Industry Support and Guidance Plan of Colleges in Gansu (2020C-20), and the Article Processing Charge (APC) was funded by Outstanding Young Scientists Support Program of Lanzhou University of Technology (LUT).

**Conflicts of Interest:** The authors declare no potential conflicts of interests with respect to the research, authorship, and/or publication of this article.

#### **References**


© 2020 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/).

## *Article* **Optimization Design of a Two-Vane Pump for Wastewater Treatment Using Machine-Learning-Based Surrogate Modeling**

#### **Sang-Bum Ma <sup>1</sup> , Sung Kim <sup>1</sup> and Jin-Hyuk Kim 1,2,\***


Received: 28 July 2020; Accepted: 15 September 2020; Published: 17 September 2020

**Abstract:** This paper deals with three-objective optimization, using machine-learning-based surrogate modeling to improve the hydraulic performances of a two-vane pump for wastewater treatment. For analyzing the internal flow field in the pump, steady Reynolds-averaged Navier-Stokes equations were solved with the shear stress transport turbulence model as a turbulence closure model. The radial basis neural network model, which is an artificial neural network, was used as the surrogate model and trained to improve prediction accuracy. Three design variables related to the geometry of blade and volute were selected to optimize concurrently the objective functions with the total head and efficiency of the pump and size of the waste solids. The optimization results obtained by using the model showed highly accurate prediction values, and compared with the reference design, the optimum design provided improved hydraulic performances.

**Keywords:** two-vane pump; Computational Fluid Dynamics (CFD); Reynolds-averaged Navier-Stokes (RANS); optimization; machine learning

#### **1. Introduction**

Recently, with the increase in the usage of disposable masks because of the COVID (Corona virus disease)-19 pandemic, used masks are being commonly discarded in toilets. If the cloth wastes such as a disposable mask or large waste such as a baby diaper is disposed into a toilet, the flow path of the pump that transports wastewater is blocked, and as a result, the function of the wastewater transport system is lost. Therefore, the demand for special pumps for wastewater transportation is increasing, and it has gained attention as an industry with the potential for future growth.

As an example of the special pumps, grinder and vortex pumps are widely used for transporting wastewater. However, these special pumps have low efficiency and high maintenance costs, contributing to large operating costs. Several studies explored the treatment of sewage containing solid waste to solve these problems [1–4]. Lu et al. [1] studied the hydraulic performance and pressure fluctuation characteristics of a grinder pump. Through a numerical analysis, the hydraulic performances of the pump when the flow path is in a clogging state and in a normal operating state were compared, and steady and unsteady Reynolds-averaged Navier-Stokes (RANS) analyses was performed. It was found that, as the degree of clogging of the grinder cutter increases, the total head of the pump declines parabolically, with the best efficiency point shifting to the low flow rate region and the high efficiency area narrowing. Schivley and Dussourd analyzed and designed a vortex pump, using a one-dimensional analytical model [2]. They compared the calculated performance parameters with

those measured by using the laboratory model, and they computed the overall hydraulic characteristics of the pump and compared these characteristics with those of many test pumps. They improved the theoretical formula used in the preliminary design of a vortex pump. Ohba et al. [3] theoretically analyzed the flow characteristics inside a vortex pump and secured the theoretical reliability by comparing the predicted values with experimental results. Their theoretical formula could predict not only the hydraulic performance of the vortex pump, but also the velocity component inside the pump. The obtained results were in good agreement with the experimentally measured values. Steinmann et al. [4] analyzed the internal flow of a vortex pump through numerical analysis and experiments to investigate the unsteady cavitating flow of a vortex pump. The Rayleigh–Plesset cavitation model was used to simulate cavitation under the overload condition of the vortex pump, and an acrylic glass window was installed in the experimental apparatus, to observe this phenomenon. Under the best efficiency point condition, the numerically derived head and shaft power of the pump were about 6% higher than the experimentally measured values. Conversely, cavitation under overload was observed more in experimental results than in numerical analysis.

In addition to the grinder and vortex pumps mentioned above, single-channel pumps designed to transport relatively larger solid wastes were actively studied [5–7]. The single-channel impeller with a single free annulus passage can smoothly transfer sewage-containing solid wastes. However, this impeller has an unsymmetrical annulus flow passage, and it is difficult to stabilize the unsteady flow-induced vibration, due to the interaction between the rotating impeller and stationary volute [5]. Shi and Tsukamoto [6] numerically analyzed the pressure fluctuation due to the impeller–diffuser interaction in a diffuser pump. They confirmed that the flow characteristics due to this interaction can be analyzed through an unsteady flow analysis. Feng et al. [7] analyzed the unsteady flow characteristics between the impeller and diffuser of a radial pump by unsteady RANS (URANS) analysis and laser Doppler velocimetry (LDV); they identified two types of rotor–stator interaction effects. One is the downstream effect induced by the impeller, which has an unsteady flow characteristic because of the highly distorted flow field and the wake of the impeller. The other is the upstream effect induced by the stator, which causes unsteady pressure and velocity fluctuations. Such a single-channel pump has the advantage of being able to transport relatively large waste solids, but it has the disadvantage that the fluid-induced vibration is greater than that in the existing special pumps (e.g., grinder and vortex pumps) because of its asymmetric structural characteristics.

The special pumps for wastewater treatment introduced so far clearly have advantages and disadvantages, depending on their type. Grinder and vortex pumps have low vibration during operation, but the size of transportable solid matters is relatively small, and the maintenance costs are relatively high. On the other hand, a single-channel pump can easily transport large solid matters, but in some cases, relatively severe vibration occurs. To solve these problems, the authors intended to design a two-vane pump in this study. The impeller of the two-vane pump is composed of two blades that are symmetrical in the rotational axis. Therefore, it is structurally simple compared to the grinder and vortex pumps, and the relatively large flow path of the impeller has a small number of blades, so waste solids can be transferred smoothly. In addition, due to the symmetrical impeller geometry, the fluid-induced vibration is relatively less than that of a single-channel pump.

In the past, design optimization using numerical analysis has been widely used for fluid-based machinery [8–11]. For example, Lee et al. [10] performed an optimization to improve the efficiency of a low-speed axial flow fan, using a gradient-based search algorithm. Lee and Kim [11] optimized axial compressor blades to improve the efficiency, using numerical optimization techniques such as conjugate direction methods and the golden section method, combined with a three-dimensional (3D) thin-layer Navier-Stokes solver. Recently, with the rapid advances in computing resources, optimal design is being actively researched based on machine-learning techniques [12,13].

In this study, a two-vane pump was designed to develop a series of special pumps for transporting wastewater. Compared to a single-channel pump, this pump can transport relatively smaller waste solids but experiences less vibration during operation because of the axial symmetry of

the impeller. Considering the characteristics of this pump, a three-objective optimization design based on machine learning was performed to maximize the size of the waste solids that can be transported while simultaneously improving the hydraulic performances of the pump. For the three-objective optimization, the geometric design variables, i.e., inlet and outlet blade angles and cross-sectional area of the volute, were chosen, and the volume of the waste solids and the head and efficiency of the pump were considered as objective functions. Their relationship was predicted by using an artificial neural network (ANN) [14] based on machine learning. Then, the genetic algorithm (GA) [15] was used to find the optimal solutions, and the Pareto-optimal front surface [16] was derived as the final optimization result.

#### **2. Numerical Methods**

#### *2.1. Two-Vane Pump Model*

The preliminary two-vane pump model used in this work was designed by using CFturbo [17], under the following design conditions: a flow rate of 0.5 m<sup>3</sup> /min, a total head of 10 m, and a rotational speed of 1760 rpm (revolution per minute).

The two-vane pump used in this study was designed for commercialization, and the preliminary design was carried out to satisfy the "KS B 6301 Standard" that requires the performance certification on fresh water under 35 ◦C, as the national certification system for the centrifugal pumps, including sewage types in the Republic of Korea. Therefore, the preliminary design and optimization of the pump were conducted to satisfy the design specifications when the working fluid is fresh water.

In this study, the initial impeller of the two-vane pump had two blades, as shown in Figure 1, and a diameter of 207 mm. The blades are shrouded impeller blades, and there is no tip clearance. The inlet and outlet blade angles are 6.50◦ and 8.00◦ , respectively, and have identical values from the hub to the shroud. To prevent the special case that the textile or cloth material caught on the blade, the ellipse ratio of the leading and trailing edge of the blade is designed to 1.0. Figure 1b shows the volute geometry and design constraints. Since the pump used in this study was installed through a manhole of the city water and sewage system, there are some constraints on the overall size of the pump and outlet diameter, as shown in Figure 1b.

The preliminary model has a specific speed (*Ns*) of 221.3 in SI (International system of units) units (rpm, m<sup>3</sup> /min, m) under the design conditions. The specific speed formula used in this study is as follows:

$$N\_s = \frac{N \times \sqrt{Q}}{\left[H\_l\right]^{3/4}}\tag{1}$$

where *N*, *Q*, and *H<sup>t</sup>* denote the rotational speed, flow rate, and total head of the pump, respectively. Figure 2 is a diagram for helping designer to roughly judge the theoretical efficiency according to the specific speed of the pump in SI unit. Now, the type of the pump can be determined by the specific speed, as shown in Figure 2, and the efficiency can be estimated by using the following approximation formula [18]:

$$\eta = \left[ 0.94 - 0.294264 \times \left[ \frac{Q}{N} \times X \right]^{-0.21333} - 12.893 \times \left[ \log\_{10} \left( \frac{2286}{N\_{\rm s}} \right) \right]^2 \right] \times 100 \text{ (\%)}\tag{2}$$

$$\mathcal{X} = \left[\frac{3.56}{\varepsilon}\right]^2\tag{3}$$

where ε is the absolute roughness height of the pump surfaces. It depends on the production processes and materials. In this study, the milling and die-casting methods that are commonly used in pump processing were assumed, and, accordingly, ε was calculated as 7.84 µm [18]. That is, according to Equations (2) and (3), the efficiency under design condition of the preliminary model is predicted to be about 71%.

(**a**) Three-dimensional (3D) geometry of the two-vane pump.

(**b**) Design constraints of the volute.

**Figure 1.** Geometry of the two-vane pump.

ଶ **Figure 2.** Efficiency according to the specific speed for the pump.

X = 3.56 ൨ *ε ε* μ ω ω ε However, based on the numerical analysis of the preliminary model, the total head was predicted as about 18 m at the design flow rate, and the specific speed was calculated as 142; these values do not satisfy the design specifications. In addition, for a specific speed of 142, the efficiency will be about 65%, as shown in Figure 2. To solve this problem, by applying the similarity law, the impeller diameter was reduced by about 13%, compared to the preliminary model, and optimization was performed to improve the hydraulic performances.

#### *2.2. Numerical Analysis*

The 3D RANS equations were solved by using a k-ω-based shear stress transport (SST) turbulence model [19] for the hydraulic analysis of the two-vane pump. The SST turbulence model is known to be suitable for predicting flow separation due to an adverse pressure gradient near the wall, and the k-ω turbulence model [20] and the k-ε turbulence model [21] are applied near the wall and to the bulk

<sup>−</sup> ω ω

flow region, respectively. These two turbulence models are connected by the blending function that is influenced by the y<sup>+</sup> value—the dimensionless number representing the distance between the wall and the first node of the grid system [19].

Commercial code ANSYS CFX 19.1 [22] was used for the RANS analysis. The computational domain is shown in Figure 1a. The 3D geometry of the impeller was created by using ANSYS Blade-Gen [22], and the volute was created by using Solidworks 2016 [23]. ANSYS TurboGrid and ICEM [22] were used to generate the computational grids for the rotating and stationary domains, respectively. The stage (or mixing plane) method was applied at the interface between the rotating and stationary domains to calculate a steady-state solution for the problems of multiple reference frames [22].

The working fluid was water at 25 ◦C. The total pressure was set to 1 atm as the atmosphere condition at the inlet boundary. The numerical analysis was performed by assigning the mass flow rate to the outlet boundary. The blade and volute surfaces in the computational domain were considered to be hydraulically smooth under an adiabatic and no-slip condition.

The grid system used in the present study consists of hexahedral grids in the rotating domain and tetrahedral grids in the stationary domain, as shown in Figure 3. In the rotating domain, O-type grids were constructed around the blades. To determine the convergence of the numerical calculations, the root-mean-square values of the residuals of the governing equations were set to be less than 10−<sup>5</sup> . The physical time scale was set to 1/ω, where ω is the angular velocity of the impeller. The computation for the steady RANS analysis was performed, using an Intel Xeon CPU with a clock speed of 2.70 GHz, and converged solutions were obtained after 1000 iterations with a computational time of approximately 8 h.

**Figure 3.** Grid systems of the two-vane pump **Figure 3.** Grid systems of the two-vane pump.

#### **3. Optimization Techniques**

The three-objective optimization problem was defined as follows:

≤ ≤ ∈ Maximize: **F**(**x**) = [F1(**x**), F2(**x**), F3(**x**)] Design variable bound: **LB** ≤ **x** ≤ **UB**, **x**∈**R**,

where **F**(**x**) is the vector of real-valued objective functions; **x** is the vector of the design variables; **LB** and **UB** denote the vectors of the lower and upper bounds, respectively; and **R** is a real number [24].

Figure 4 shows the procedure of a typical optimization design. Once the operating conditions of the design target are determined, the type of turbomachinery and airfoil (or hydrofoil) are selected. Then a preliminary design is performed through one-dimensional mean-line analysis, and the initial blade topology is derived. Subsequently, the optimization design is performed for improving the performance of the turbomachinery.

**Figure 4.** Procedure of the optimization design based on machine learning.

In this study, the optimization design was based on machine learning. First, the objective functions and constraints were defined according to the design goal. Subsequently, the design variables and their ranges were chosen. Thereafter, a database of the correlation between the design variables and the objective functions was established within the design space, using the design of experiment (DOE). In the next step, a predictive model was constructed by using ANN [14] to correlate the design variables and objective functions of the two-vane pump. In this step, the predictive model was trained by using machine learning. This process is described in detail in Section 3.2. The GA [15] was used to find the optimal design solution, considering the correlation of each objective function in the constructed prediction model. The GA is a well-known stochastic searching algorithm based on the mechanism of natural selection, genetics, and evolution. It evaluates various points in the design space and can be applied to find a global solution to any given problem. This algorithm proceeds as shown in Figure 5. Finally, the Pareto-optimal front surface was derived by using MATLAB R2018b [25], and the optimization procedure was terminated.

*η*

*η*

**Figure 5.** Genetic algorithm. **Figure 5.** Genetic algorithm.

*ρ τ ω*

*ρ τ ω*

(௨௧௧ − ௧) 

> <sup>௦</sup> = 1 6 <sup>௦</sup> ଷ

௨௧௧ − ௧ 

× 100%

<sup>௧</sup> =

=

#### *3.1. Design Variables and Objective Functions*

In order to maximize the hydraulic performances and size of waste solids, the efficiency (η) and total head (*Ht*) of the pump and the volume of waste solids (*Vs*) were selected as the objective functions:

$$H\_t = \frac{P\_{outlet} - P\_{inlet}}{\rho \text{g}} \tag{4}$$

$$\eta = \frac{(P\_{outlet} - P\_{inlet})Q}{\tau \omega} \times 100\% \tag{5}$$

$$V\_s = \frac{1}{6}\pi D\_s^3\tag{6}$$

where *P*, ρ, *g*, *Q*, τ, and ω are the total pressure, density of the working fluid, acceleration due to gravity, flow rate, torque, and angular velocity, respectively. Further, *D<sup>s</sup>* in Equation (6) is defined as the distance between the leading edge of the blade and the other blade at the impeller inlet, as shown in Figure 6.

**Figure 6.** Definition of the waste solid volume.

*β β* Figure 7 shows the design variables considered in this work. Their ranges are listed in Table 1. The inlet and outlet blade angles and the cross-sectional area of the volute were chosen for the optimization. They are defined as shown in Figure 7. As mentioned earlier, the blade angle distribution is the same in the span direction from the hub to the shroud and is defined by using the fourth-order Bézier curve [26] in the streamwise direction, as shown in Figure 7a. In order to adjust the inlet (β1) and outlet (β2) blade angles, two control points (CP<sup>1</sup> and CP2) were fixed. The distribution of the cross-sectional area of the volute was changed linearly in the circumferential direction, as shown in Figure 7b. *β β*

(**c**) Definition of theta.

**Figure 7.** Definition of the design variables.


**Table 1.** Ranges of design variables.

LB = lower bound; UB, upper bound.

#### *3.2. Surrogate Modeling Based on Machine Learning*

In the present optimization design, supervised machine learning was adopted, considering the reasonable computation costs. This optimization technique is a very efficient approach for optimizing a system without analytical representation, fitting a surrogate model.

In order to construct the input data, the central composite design was used as the DOE. In all, 15 experimental points were arranged for three design variables, and the objective function values were derived through RANS analysis at each design point. These values are listed in Table 2. In this table, the values of the design variables and objective functions are normalized by dividing them by the corresponding reference value.


**Table 2.** Initial input data for the supervised machine learning.

A radial basis neural network (RBNN) model [27], which is a type of ANN, was used in this optimization study. The RBNN model has a hidden layer of radial neurons and an output layer of linear neurons, as shown in Figure 8. The hidden layer uses a series of radial primitives to nonlinearly modify the input space to the intermediate space. The output of the hidden layers then executes a linear combiner to produce the desired targets [27]:

$$f(\mathbf{x}) = \sum\_{j=1}^{N} w\_j \phi\_j \tag{7}$$

where *w<sup>j</sup>* is the weight, and ϕ*<sup>j</sup>* is radial basis function, which uses the Gaussian function. Several parameters are needed to construct a surrogate model: the input weight and the center and width of a unit for the Gaussian function. In the present RBNN model, the input weights were chosen by machine learning. The input vector with the worst performance was chosen as the center of a new hidden-layer Gaussian function [25]. Then the RBNN only needed to determine the width of the Gaussian function (spread constant). The network training was performed by adjusting the cross-validation error (CV) by changing the spread constant (SC), as shown in Figure 9. SC1, SC2,

and SC<sup>3</sup> correspond to the total head, efficiency, and size of waste solids, respectively. The SC values for each objective function were chosen by a K-fold CV test [28].

(x) =

*φ*

ே

ୀଵ

*β β β β*

**Figure 8.** Schematic of the radial basis neural network.

(**c**) Spread constant for Vsv.

**Figure 9.** Cross-validation errors vs. spread constant (SC) values.

(( , ); = 1, … , N) The data sample ((*x<sup>j</sup>* , *yj*); *j* = 1, . . . , N) was partitioned into K disjoint subsets (K-fold CV), as shown in Figure 10. Of these, (K-1) folds were used to train the RBNN network, and the last fold (the Kth set) was used for evaluation. This process was repeated K times, using a different fold for evaluation each time. The network training was performed by adjusting the CV error by changing SC. The CV error at a particular SC value was calculated as follows:

*β /β β /β / η η V /V*

ଶ

ୀ௧ ௦௧ ,

CV (SC) =

ଵ 

∈ ො

∑ ∈ 

*Processes* **2020**, *8*, 1170

$$\text{CV (SC)} = \frac{1}{\text{K}} \sum\_{i=1}^{k} \in\_{k\prime} \in\_{k} = \sum\_{j=k\text{th set}} (y\_j - \hat{y}\_j)^2 \,. \tag{8}$$

where ∈*<sup>k</sup>* is the prediction error for the Kth set. The predicted values *y*ˆ*<sup>j</sup>* were determined by using the constructed RBNN model from the sample points in the (K-1) subsets. In the present study, K was set as 15. According to the results of the K-fold CV test, the final SC values SC1, SC2, and SC<sup>3</sup> were set as 0.3, 1.9, and 9.9, respectively.

**Figure 10.** K-fold cross-validation.

#### **4. Results**

#### *4.1. Grid Dependency Test and Validation of Numerical Results*

To eliminate the grid dependency of the numerical solutions, grid dependency tests were performed in the computational domain, as shown in Figure 11. In these tests, the efficiency and total head under the design condition were compared, and their values were normalized by dividing by the corresponding convergence values. As shown in Figure 11, the grid system with 3.4 <sup>×</sup> <sup>10</sup><sup>6</sup> nodes undergoes only 0.002% and 0.083% changes in the efficiency and total head, respectively, compared with the grid system with 3.0 <sup>×</sup> <sup>10</sup><sup>6</sup> nodes. Hence, the latter grid system was selected as the optimal grid system for the computational domain.

**Figure 11.** Grid dependency test.

In earlier research, the authors analyzed the internal flow characteristics of hydraulic machines, such as pumps and water turbines, through numerical analysis and compared and verified the numerical results through experiments [29–32]. In a previous study, a single-channel pump with design specifications similar to those of the two-vane pump considered in this study was developed, and Figure 12 shows the results of the experiments conducted in the previous study [29]. The hydraulic performances, i.e., the head coefficient and efficiency, were compared with the corresponding experimental data in the operating ranges. The total head values derived from the numerical analysis were almost identical to the experimental data. Meanwhile, the efficiency values were relatively higher than those of the experimental data because the numerical analysis did not include the mechanical losses. However, the general trend of the efficiency curve matched well. In this study, the numerical scheme used in the previous studies of the authors [29–32] was adopted, and the numerical analysis technique was verified through several peer reviews. The results derived of this study will be verified through experiments in future work.

(**a**) Prototype model for experiments. (**b**)Verification results of the prototype model.

**Figure 12.** Validation data of previous study [29].

#### *4.2. Optimization Results*

Figure 13 shows the Pareto-optimal front surface for the three-objective optimization. This surface was obtained by using the GA based on the RBNN predictive model. To investigate the accuracy of the optimization results, five arbitrary optimal designs (AODs) were selected. They are listed in Table 3 and shown in Figure 13.

**Figure 13.** Pareto-optimal front surface with arbitrary optimal designs (AODs).


**Table 3.** Optimization results.

RANS = Reynolds-averaged Navier-Stokes.

The objective functions of the predicted AODs were calculated by RANS analysis and compared with the predicted objective function values listed in Table 3. The RANS results indicate that the maximum relative errors were less than 1.90% for the total head, 1.80% for the efficiency, and 4.00% for the size of waste solids. Thus, the surrogate model is regarded as being constructed based on reliable data, and the results of optimization indicate the excellent accuracy of surrogate prediction.

All the AODs are predicted to have a normalized total head of more than 0.611, which satisfies the design target, and the efficiency and the size of waste solids are inversely related. In addition, all the AODs have less cross-sectional area of the volute, compared to the reference design, and the efficiencies of the AODs are improved. The design focused on the size of waste solids (AOD 1) doubles the solid size compared with that in the reference design, whereas the design focused on efficiency (AOD 5) shows an increment of 14.14% in the efficiency.

#### **5. Discussion**

The size of waste solids is directly related to the blade inlet angle. As the blade inlet angle increases, the size of waste solids increases. The flow passage areas of the AODs in the meridional direction are compared in Figure 14. Their values are normalized by using the maximum area of AOD 1 in the meridional direction. It was confirmed that the inlet passage area of AOD 1 with the greatest increase in the blade inlet angle was the largest, and the inlet passage area of AOD 5 with the reduced inlet blade angle was the smallest among AODs.

**Figure 14.** Flow passage area distribution in the meridional direction.

AOD 3, which satisfies the design target, was selected for further analysis of the optimum design. The 3D geometry for AOD 3 is shown in Figure 15. Through optimization, an inlet blade angle of about 1.289 times the reference design and an outlet blade angle of about 0.128 times the reference design were realized. That is, considering the change in the flow path inside the impeller, the blade wrap angle was reduced compared to the reference design. Reduction in the blade wrap angle led to a shorter flow path, thus reducing the friction loss. In the case of the volute, the cross-sectional area of the optimum design was less than that in the reference design. Furthermore, the cross-sectional area continuously decreased from the tongue to the outlet, and then it decreased by approximately 18.9% at the outlet, as shown in Figure 16.

**Figure 15.** Comparison of three-dimensional (3D) geometries.

**Figure 16.** Cross-sectional area distribution of the volute.

Figure 17 shows the velocity distribution at 50% span of the blade. In the optimum design, the stagnation point is formed at the leading edge of the blade, whereas in the reference design, it is located relatively downstream. This phenomenon occurs because the inlet angle of the blade is not designed to fit the operating condition, and it can increase the incidence angle to cause flow separation and reduce the hydraulic performance of the pump. In addition, the stagnation point is formed on the suction side (SS) of the blade, and the flow proceeds non-uniformly. In particular, a very high-velocity region is distributed at the leading edge (LE) of the pressure side (PS), and severe flow separation occurs at the SS, resulting in a blockage inside the passage.

**Figure 17.** Velocity distributions at 50% span.

*β β* *β*

The distribution of the streamlines and the inlet velocity component at 50% span are analyzed and shown in Figure 18. As described in Figure 17, the inlet flow collides at the SS of the blade in the case of the reference design; conversely, in the optimum design, the inlet flow collides precisely with the LE of the blade. Through the optimization, the blade inlet angle (β*<sup>b</sup>* ) of the optimum design is greater than that in the reference design, and, therefore, the inlet flow proceeds smoothly, as shown in Figure 17b. *β*

(**b**) Optimum design.

**Figure 18.** Streamlines at 50% span.

−

*β β* The rotational velocity (*U*) of the blade, the axial velocity (*Ca*), and the relative velocity (*W*) of the flow at the inlet of the impeller can be represented in the velocity triangle diagram. When the three velocity components are known, the flow angle (β) can be derived, and the incidence angle (*I*) can be calculated by comparison with the blade angle (β*<sup>b</sup>* ). The incidence angle distribution from the hub to the tip is shown in Figure 19. The incidence angle of the reference design gradually decreases from the hub to about 80% span and then increases again to the tip span. The tendency of the optimum model is similar to that of the reference design, but the incidence angle is small overall in all spans, compared with the reference design. The largest incidence angle in the reference design is 24◦ at the tip region, whereas that in the optimum design is 18◦ at 15% span.

**Figure 19.** Incidence angle distribution at the leading edge (LE) along the span.

In rotating-fluid machines, the flow is driven by the centrifugal force. Hence, the incidence angle at the tip region is smaller than that at the hub; this improves the hydraulic performance. In particular, since the pump used in this study is a special pump for transporting wastewater, the blade angle distribution from the hub to the tip should be maintained. From this point of view, it is judged that the incidence angle distribution of the optimum design shown in Figure 19 is ideal.

Figure 20 shows the distribution of streamlines and vortices on the iso-surface of the velocity invariant (><sup>5</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> s −2 ), to visualize the flow characteristics inside the impeller. As shown in Figures 17 and 18, in the case of the reference design, the inlet flow angle is significantly different from the blade inlet angle, and an apparent inter-blade vortex is observed inside the passage, as shown in Figure 20a. On the other hand, the inlet flow angle of the optimum design agrees well with the blade inlet angle, and the flow proceeds smoothly into the passage, as shown in Figure 20b.

− **Figure 20.** Three-dimensional (3D) streamlines with vortices distribution on the iso-surface of the velocity invariant (><sup>5</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> s −2 ).

In order to analyze the blade loading, the pressure distribution on the blade in the optimum and reference design are compared, as shown in Figure 21. The pressure values are normalized by using the maximum pressure value of the reference design. Similar to Figures 17 and 18, Figure 21 shows that the static pressure of the reference design at LE is relatively less than that in the optimal design, owing to the rapid increase in the velocity in the PS of the blade LE. In addition, the pressure on the PS and SS of the reference design shows a larger overall region than in the optimum design; in other words, the blade loading of the reference design is larger.

−

**Figure 21.** Pressure distributions at 50% span.

Figure 22 shows the internal velocity distribution from the volute tongue to the outlet. In the reference design, a relatively high-velocity region is distributed near the tongue, but the velocity decreases toward the outlet. On the other hand, the optimum design shows a relatively uniform velocity distribution over the entire area of the volute. In the reference design, the cross-sectional area of the volute is designed to be larger than necessary, and as the flow diffuses near the outlet, it is judged that the static pressure recovery occurs more than in the optimum design.

**Figure 22.** Velocity distributions inside the volute.

Figure 23 shows the iso-surface contours for a low velocity of 1.0 m/s, to identify the low-velocity region inside the pump. As shown in Figure 22, a relatively wide low-velocity area is formed at the outlet of the reference design. Since the two-vane pump considered in this work is used for transporting wastewater, the solids move inside the impeller and volute. Now, the presence of low-velocity region can cause solids to stagnate. Therefore, the pump should be designed so that no such low-velocity region is created. From this point of view, the optimum design is an ideal design with a significant reduction in the area of the low-velocity region.

**Figure 23.** Iso-surface for the low-velocity region (<1 m/s).

#### **6. Conclusions**

In the present study, the shapes of a blade and volute were optimized to improve the hydraulic performance of a two-vane pump for transporting wastewater via 3D RANS analysis based on supervised machine-learning optimization. Considering the waste solids, the blade inlet and outlet angles and the cross-sectional area of the volute were selected as the design variables. Three objective functions were considered to maximize the hydraulic performances, i.e., the total head and efficiency of the pump, while maximizing the size of the waste solids.

In this analysis, the RBNN predictive model was trained by using machine learning, and the GA was used to find global optimal solutions. The network training was performed by adjusting the CV error by changing the SC. Through training, an accurate and reliable RBNN surrogate model was constructed, and five arbitrary optimum design points were derived. They were compared with the numerical results to verify the predicted accuracy. The machine-learned surrogate model showed very accurate predictions compared with the numerical results and had relative errors less than 5%.

Through the optimization, the pump was made compact, and the efficiency was improved by about 14% compared to the reference design, while satisfying the design goal of the total head (>10 m). As the cross-sectional area of the optimum volute was reduced, excessive diffusion at the outlet was reduced, and, therefore, the velocity distribution inside the volute was more uniform compared with the reference design. In addition, as the inlet angle of the optimum blade increased, the water flowing from the impeller inlet proceeded smoothly into the flow path, and the flow separation and formation of the inter-blade vortex were significantly reduced.

In future work, the hydraulic performances of the optimum design will be investigated in more detail experimentally, and the characteristics of the flow, containing solid wastes, will be also analyzed.

**Author Contributions:** J.-H.K. and S.K. designed the two-vane pump; S.-B.M. and S.K. analyzed the numerical data; S.-B.M. contributed the optimization techniques and analyzed optimization results and organized draft paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Demand-Based-Manufacturing Technique Commercialization R&D Project of the Korea Institute of Industrial Technology (KITECH) (No. JB200007), which was funded by the Ministry of Science and ICT (MSIT). The authors gratefully acknowledge this support.

**Conflicts of Interest:** The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

#### **References**


© 2020 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/).

## *Article* **Applications of an Improved Aerodynamic Optimization Method on a Low Reynolds Number Cascade**

### **Shuyi Zhang, Bo Yang \*, Hong Xie and Moru Song**

School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China; shuyizhang@sjtu.edu.cn (S.Z.); xiehong1211@sjtu.edu.cn (H.X.); samuel0624@sjtu.edu.cn (M.S.) **\*** Correspondence: byang0626@sjtu.edu.cn

Received: 18 August 2020; Accepted: 9 September 2020; Published: 14 September 2020

**Abstract:** The effect of cascade aerodynamic optimization on turbomachinery design is very significant. However, for most traditional cascade optimization methods, aerodynamic parameters are considered as boundary conditions and rarely directly used as the optimization variables to realize optimization. Given this problem, this paper proposes an improved cascade aerodynamic optimization method in which an incidence angle and nine geometric parameters are used to parameterize the cascade and one modified optimization algorithm is adopted to find the cascade with the optimal aerodynamic performance. The improved parameterization approach is based on the Non-Uniform Rational B-Splines (NURBS) method, the camber line superposing thickness distribution molding (CLSTDM) method, and the plane cascade design method. To rapidly and effectively find the cascade with the largest average lift-drag ratio within a certain range of incidence angles, modified particle swarm optimization combined with the modified very fast simulated annealing algorithm (PSO-MVFSA) is adopted. To verify the feasibility of the method, a cascade with NACA4412 and a practical cascade are optimized. It is found that the average lift-drag ratios of two optimal performance cascades are respectively increased by 13.38% and 15.21% in comparison to those of two original cascades. Meanwhile, through optimizing the practical cascade of the Blade D500, under different volume flow rates, the pressure coefficient of the optimized cascade is increased by an average of more than 6.12% compared to that of the prototype, and the average efficiency is increased by 11.15%. Therefore, this improved aerodynamic optimization method is reliable and feasible for the performance improvement of cascades with a low Reynolds number.

**Keywords:** cascade; aerodynamic; parameterization; plane cascade design; incidence angle; PSO-MVFSA; optimization

### **1. Introduction**

Blade design is of great importance to the efficiency and properties of turbomachinery. Achieving the aerodynamic design of a blade is a very complex and arduous task due to complicated flow phenomena and the interactions among various parameters [1]. A small geometric change of one blade can lead to a deterioration of the aerodynamic performance of the whole machine [2]. In general, the cascade design method is one of the most popular design methods employed for turbomachinery [3,4]. In this method, the blade is stacked by the sections on the different radiuses and shown in Figure 1a. Additionally, the section is projected onto the plane to form the cascade, as shown in Figure 1b. Therefore, the blade performance is affected by the cascade design. With the rapid development of the computer technique, more and more parameterization methods and optimization algorithms have been proposed and used in the design process of a cascade. This means that the time required to design a new cascade is becoming shorter and the aerodynamic performance of the new cascade is becoming better.

**Figure 1.** An explanation of blades and a cascade.

A good parameterization method can not only use fewer design variables to describe an airfoil accurately, but also rapidly re-construct one airfoil in the optimization process [5,6]. The parameterization methods are usually divided into two categories: The constructive method and the deformative method. Each method has been continuously improved by many researchers [7]. The deformative parameterization method is the simpler of the two methods. In this method, a standard airfoil is deformed to generate one new airfoil, in order to satisfy a certain condition. The Hicks-Henne function [8], radial basis function [9], Bezier function [10], B-Spline function [11], and Non-Uniform Rational B-Splines (NURBS) function [12–14] are usually used as deformative functions to generate a new airfoil based on the original airfoil. In particular, NURBS [12] is the most popular function due to its ability of local control and its conics description over the curve. However, these deformative functions are used to generate one new airfoil based on the point coordinates of an airfoil. To relate the airfoil shape to the airfoil geometric feature parameters, the camber line superposing thickness distribution molding (CLSTDM) method [15] was proposed. In this method, several airfoil geometric parameters are used to parameterize the half-thickness distribution curve and the mean camber curve through two polynomials. Then, these two curves are coupled to form a whole airfoil. The feasibility of this method has been proved by several works [15–17]. In the CLSTDM method [15,16], the blade contours described by many coordinate points can be transformed into functions controlled by several parameterized variables. Moreover, it is convenient for designers to use this method to parameterize one blade based on their experience. Nowadays, it is widely used in the optimization of turbomachinery. However, the aerodynamic parameters are not considered in the parameterization and still act as the boundary condition of flow field computing.

on

To rapidly find the airfoil with the optimal aerodynamic property, many intelligent optimization technologies have been used in the process. Among all of the published optimization algorithms, the genetic algorithm (GA) [18–21] and the simulated annealing (SA) algorithm [10,11,22], as two traditional intelligent optimization algorithms, have been widely used in airfoil optimization. They aim to find the airfoil with the optimal performance precisely. However, much time is required to complete the search, which leads to a poor computational efficiency [23,24]. A new intelligent algorithm, known as particle swarm optimization (PSO) and proposed by Kenney and Eberhart [25], can be used to solve this problem. PSO is a population-based, self-adaptive searching optimization method. The principle of this method is based on animal social behaviors, such as birds' migration. However, it is easy to become trapped into a local extreme value or converged to precocity by the standard PSO. Therefore, researchers began looking for improved methods to solve this problem. Shi [26] used a linearly decreasing inertia weight to balance the global and local searching character. However, the local searching capability of this method was weak. Simultaneously, it is hard to predict the maximum iteration number. Clerc [27] set a constriction factor determined by two learning factors to cancel the boundary limits of velocity, and to balance the global and local searching capability. Hu [28] adopted a stochastic inertia weight to replace the linearly decreasing inertia weight. This method could accelerate the convergence velocity

and avoid being trapped into a local best solution. However, these methods still have the risk of local convergence.

Most approaches to parameterization of the cascade have been used with only the help of geometric feature parameters, and the aerodynamic parameters were not referred to. The efficiency of optimization was thus negatively affected. Therefore, in this paper, considering the cascade aerodynamic characters, one aerodynamic parameterization approach for a low Reynolds number cascade is proposed. To find the airfoil with the optimal performance during a certain range of incidence angles, a modified PSO-MVFSA algorithm is studied. Furthermore, two cases, such as the cascade with NACA4412 and the blade of FAN D500, are selected to verify the feasibility of the improved parameterization and optimization method.

#### **2. Aerodynamic Parameterization Method**

#### *2.1. CLSTDM-NURBS Method*

=1.0

,

In this paper, for the CLSTDM method, the pressure side and suction side are obtained through the camber curve superposing the half-thickness distribution curve. This method is used to describe an airfoil and is shown in Figure 2. Due to its good ability to design a complex geometry, the two-order NURBS function defined by Equation (1) is used to describe the camber curve and the half-thickness distribution curve.


$$\begin{cases} \mathcal{C}(u) = \sum\_{i=0}^{n} \frac{B\_{2,i}(u)\omega\_i Q\_i}{B\_{2,i}\omega\_i(u)}\\\ B\_{2,i}(u) = \frac{2!}{i!(2-i)!}u^i(1-u)^{2-i} \end{cases} \tag{1}$$

χ2 -

2, 2 2, 2! () 1 !2 ! , =2 0,1,2 ( ) where *u* is the knot vector, *n* is the order of NURBS (*n* = 2), *i* is the mark of the control points (*i* = 0, 1, 2), *C*(*u*) is the coordinate of the point of the fitting curve parameterized by the NURBS function, *B*2,*i*(*u*) is the Bernstain function, ω*<sup>i</sup>* represents the weight coefficients (ω0= 1, ω2= 1, ω<sup>1</sup> = ω), and *Q<sup>i</sup>* is the control point. To solve the two-order NURBS function, the De Boor algorithm [29], which provides a fast and numerically stable way of finding a point on a B-spline curve with the given *u* in the domain, is adopted and programed.

2, ( ) 0 2 =1, =1 <sup>1</sup>= <sup>012</sup> , , 3 4 , , For camber parameterization, the camber curve is divided into two two-order NURBS curves (solid line), and these two NURBS curves are respectively controlled by several geometric control points (*Q*0, *Q*1, *Q*2, *Q*3, *Q*4) shown in Figure 3. Five geometric feature parameters of an airfoil, such as the chord line *L* = 1.0, the leading edge angle χ1, the trailing edge angle χ2, and the coordinate of the maximum camber point *Bx*, *B<sup>y</sup>* , are adopted to derive the coordinates of the geometric control points of NURBS. To ensure the continuous property of two NURBS curves at the location of the maximum camber point *Q*2, two lines *Q*1*Q*<sup>2</sup> and *Q*2*Q*<sup>3</sup> need to be collinear. Moreover, all of the geometric control points of NURBS are derived as follows in Equation (2).

$$\text{The first-order coupling between the two-dimensional } \mathcal{N} \text{-matrices is the only possible } \mathcal{N} \text{-matrices with } \mathcal{N} = \{0, 1, 2, \dots, N\} \text{ and } \mathcal{N} = \{0, 1, 2, \dots, N\}.$$

χ1 χ

1

**Figure 3.** Camber curve parameterization. 3 tan ,

2

2

Q2 (Bx , B<sup>y</sup> <sup>Q</sup> ) <sup>Q</sup><sup>3</sup>

4

$$\begin{aligned} \left[ \begin{array}{c} Q\_0\\ Q\_1\\ Q\_2\\ Q\_3\\ Q\_4 \end{array} \right] = \left[ \begin{array}{c} (0,0)\\ \left(B\_y/\tan(\chi\_1), B\_x\right) \\\left(B\_x, B\_y\right) \\\left(L - B\_y/\tan(\chi\_2), B\_y\right) \\\left(L,0\right) \end{array} \right] \\ \quad \text{(2)}\\ \text{M} \left(\text{which are similar with} \begin{array}{c} \text{minimize} \\\\ \text{ $\dot{\chi}$  and  $\dot{\chi}$  is similar.} \end{array} \right) \end{aligned} \tag{2}$$

Q<sup>4</sup> L (L , 0)

χ2

3 2 4 tan , , 0 <sup>012345678</sup> ,,,,,,,, *R*<sup>1</sup> *R*<sup>2</sup> =1.0 For parameterization of the half-thickness distribution curve, the curve is divided into three parts, including the leading edge half-thickness, the middle half-thickness, and the trailing edge half-thickness, as shown in Figure 4. The leading edge (LE) part is described by one two-order NURBS curve, the middle part by two two-order curves, and the trailing edge (TE) part by one two-order NURBS curve. Four NURBS curves are controlled by nine geometric control points (*P*0, *P*1, *P*2, *P*3, *P*4, *P*5, *P*6, *P*7, *P*8), which are derived by seven geometric feature parameters of the airfoil, such as the leading edge radius *R*1, the trailing edge radius *R*2, the chord line *L* = 1.0, the thickness gradient angles α1, α2, and the coordinate of the maximum half-thickness point (*Tx*, *Ty*). As is the case for the camber curve, the half-thickness distribution curve also requires continuity. Therefore, at the three geometric control points *P*2, *P*4, *P*6, it is necessary to maintain collinearity for the relative lines. Utilizing geometric principles, these control point coordinates are shown in Equation (3). <sup>012345678</sup> ,,,,,,,, <sup>1</sup> <sup>2</sup> =1.0 1 2 , (,) <sup>246</sup> , ,

**Figure 4.** Half-thickness distribution curve parameterization.

$$\begin{pmatrix} P\_0 \\ P\_1 \\ P\_2 \\ P\_3 \\ P\_4 \\ P\_5 \\ P\_6 \\ P\_7 \\ P\_8 \\ P\_9 \\ P\_9 \end{pmatrix} = \begin{bmatrix} (0,0) \\ (0,R\_1\*\cos(\alpha\_1)-R\_1\*\tan(\alpha\_1)\*(1-\sin(\alpha\_1)))\tan(\alpha\_1) \\ (R\_1\*(1-\sin(\alpha\_1)),R\_1\*\cos(\alpha\_1)) \\ (T\_y-R\_1\*\cos(\alpha\_1))/\tan(\alpha\_1)+R\_1\*(1-\sin(\alpha\_1)),T\_y \\ (T\_{y\_1}T\_y) \\ -(T\_y-R\_2\*\cos(\alpha\_2))/\tan(\alpha\_2)+L-R\_2\*(1-\sin(\alpha\_2)),T\_y \\ (L-R\_2\*(1-\sin(\alpha\_2)),R\_2\*\cos(\alpha\_2)) \\ (L,R\_2\*\cos(\alpha\_2)-R\_2\*\tan(\alpha\_2)\*(1-\sin(\alpha\_2))) \\ (L,0) \end{bmatrix} \tag{3}$$

, cos( ) tan( ) (1 sin( ))

cos( ) tan( ) 1 sin(

2

) ,

2 1 1 1 1

22 2 2 <sup>6</sup>

2 2 2 2

2 22 2 2

1 sin( ) , cos( ) , cos( ) tan( ) (1 sin( ))

(1 sin( )) , cos( )

2 2 2 2

3

8 9

P1

P0

4 5

= ,

, 0

, 0

7 8 9

#### *2.2. Improved Aerodynamic Parameterization*

Unlike conventional parameterization approaches in which the airfoil is only parameterized with the geometric feature parameters [30–32], an improved aerodynamic parametrization approach combining the plane cascade design method and the CLSTDM-NURBS is proposed, in which the incidence angle, *i*, and nine geometric parameters are used as control variables.

For one specific cascade, it can be assumed for the airflow condition that the inlet flow angle β<sup>1</sup> is constant and along the tangential line. In this case, a change of the incidence angle can cause a change of the geometric inlet angle of the cascade, and the variation of the incidence angle ∆*i* is equal to that of the geometric inlet angle ∆β1*A*. From the definition of the geometric inlet angle, it is clear that the slope of the mean camber curve at the leading point is changed with the changing of the geometric inlet angle. This means that the variation of the incidence angle ∆*i* has an indirect influence on the modified value of the leading edge angle χ ′ 1 . The relationship is shown in Equation (4).

$$
\chi\_1' = \chi\_1 - \Delta i \tag{4}
$$

Simultaneously, the deviation angle δ can also be affected by a change of the incidence angle. In order to determine the deviation angle, one semi-rational formula of the deviation angle was used, based on the plane cascade experiments by Howell [33,34], which was suitable for a low Reynolds number cascade and only applicable in the application conditions (τ = 0.7 ∼ 2.0, *T<sup>y</sup>* = 0.05 ∼ 0.12, and *B<sup>x</sup>* = 0.4 ∼ 0.5). The definition formulas of the incidence angle *i* and the deviation angle δ have been proposed in the cascade design process [3,35].

$$
\pi = \frac{t}{L} \tag{5}
$$

$$
\delta = \beta\_{2A} - \beta\_2 \tag{6}
$$

In this paper, one special equation combining a semi-empirical formula [33,34] is shown in Equation (7), where ϕ is equal to 0.5 for the rotor blade. Therefore, if the aerodynamic and geometric parameters *i*, τ, β1, β2, and *B<sup>x</sup>* are given, the geometric outlet angle β2*<sup>A</sup>* can be calculated by Equation (7). Additionally, the revised trailing edge angle χ ′ 2 can be obtained by Equation (8). However, in Equation (7), it is found that the cascade solidity τ and the outlet flow angle β<sup>2</sup> cannot be determined.

$$\mathbf{a} \left( \beta\_{2A} - \beta\_{2} \right) = \left[ 0.23 \left( 2B\_{\mathbf{x}} \right)^{2} - 0.002 \beta\_{2} + 0.18 \right] \left( \beta\_{2A} - i - \beta\_{1} \right) \left( \frac{1}{\tau} \right)^{\varphi} \tag{7}$$

$$
\chi\_2' = \chi\_2 - \Delta\beta\_{2A} \tag{8}
$$

In order to obtain the abovementioned two parameters, the diffusion factor *D* related to the cascade solidity τ, the inlet flow angle β1, and the outlet flow angle β2, is introduced as a constraint. This coefficient can be used to control the aerodynamic load of the cascade. The definition of the diffusion factor is shown in Equation (9).

$$D = \left(1 - \frac{\sin \beta\_1}{\sin \beta\_2}\right) + \frac{\sin \beta\_1}{2\pi} (ctg\beta\_1 - ctg\beta\_2) \tag{9}$$

Therefore, during the process of parameterization, the leading edge angle and trailing edge angle shown in Figure 3 are independent variables and no longer depend on the incidence angle. In this way, the control variables of cascade parameterization are determined, such as *Bx*, *By*, *Tx*, *Ty*, *R*1, *R*2, α1, α2, and *i*.

#### **3. The Improved Airfoil Aerodynamic Optimization Method**

#### *3.1. Modified PSO-MVFSA*

With the extensive use of standard particle swarm optimization (PSO), some drawbacks have been found, such as the local extreme minimum and the precocity. To solve these problems, many works have been published [26–28]. Some only adjusted the change of inertia weight to avoid being trapped into a local optimal solution. Moreover, some used two learning factors to balance the local searching and global searching. In this paper, the control variables of cascade parameterization are grouped as a particle. Additionally, the best particle position *Xi*,*<sup>j</sup>* corresponding to the optimal cascade is obtained by the modified PSO, as shown in Equation (10) [27].

$$\begin{cases} V\_{i,j}^1 &= \phi \cdot \left[ \omega V\_{i,j}^0 + c\_1 \gamma\_1 \left( P\_{i,j} - X\_{i,j}^0 \right) + c\_2 \gamma\_2 \left( P\_{\mathcal{G},j} - X\_{i,j}^0 \right) \right] \\ \phi &= \frac{2}{\left| 2 - (c\_1 + c\_2) - \sqrt{(c\_1 + c\_2)^2 - 4(c\_1 + c\_2)} \right|} \\ \omega &= \mu\_{\min} + (\mu\_{\max} - \mu\_{\min}) \* \gamma\_4 + \sigma \cdot \gamma\_3 \\ X\_{i,j}^1 &= X\_{i,j}^0 + V\_{i,j}^1 \end{cases} \tag{10}$$

where *V* is the particle velocity, subscript *i*, *j* is the particle sequence, ω is the stochastic inertia weight [28], φ is the constriction factor [27], *c*1, *c*<sup>2</sup> represents two learning factors, γ1, γ<sup>2</sup> represents the random number uniformly distributed in (0, 1), *Pi*,*<sup>j</sup>* is the best position in its flight history, *Pg*,*<sup>j</sup>* is the best position in the particle swarms, µmax is the maximum stochastic inertia weight, µmin is the minimum stochastic inertia weight, σ is the variance of the stochastic inertia weight, and γ3, γ<sup>4</sup> is the random number of the standard normal distribution.

However, it is very difficult to judge whether the results are the local optimal solutions or the global solutions when only using the modified PSO. Due to the probability of the simulated annealing (SA) algorithm jumping out of the local value, it can be used to solve the problem. Furthermore, considering that the low searching velocity of the standard SA algorithm can lead to a greater time consumption, modified very fast simulated annealing (MVFSA) [36] is adopted to search for the global optimal minimum. Since the Cauchy distribution depending on temperature was better than the Gaussian distribution, the coefficient of variable perturbation *yij* based on the Cauchy distribution has been redefined by Equation (11). In order to improve the efficiency further, the random probability of the acceptance, *P<sup>r</sup>* , has been redefined based on the Boltzmann–Gibbs distribution, as shown in Equation (12).

$$y\_{ij} = T\_{\text{max}} \exp\left(-ck^{1/N\_1}\right) \cdot \text{sgn}(\mu - 0.5) \left[\left(1 + \frac{1}{T}\right)^{|2\mu - 1|} - 1\right] \tag{11}$$

$$P\_r = \left[1 - (1 - h)\left(E\{M\_{ij}'\} - E\{M\_{ij}\}\right) / T\right]^{1/(1-h)},\tag{12}$$

where *u* is a random number ranging from 0 to 1, *T*max is the initial simulated high temperature, *N*<sup>1</sup> is the syllogism coefficient, *k* is the number of the marked annealing stage, *c* is a given constant value, *M*′ *ij* is the disturbed variable, *Mij* is the undisturbed variable, and *h* is a real number (set *h* = 0.5). Furthermore, *E <sup>M</sup>ij* is adopted to evaluate the current energy of the particle.

#### *3.2. Verification of Modified PSO-MVFSA*

Due to its local optimums and premature results, the Rastrigin function is usually used as a test function. The definition of the function is shown in Equation (13).

$$f(\mathbf{x}) = \sum\_{j=0}^{D} \mathbf{x}\_j^2 - 10 \cdot \cos(2\pi \mathbf{x}\_j) + 10 \qquad \qquad \qquad \mathbf{x}\_j \in [-20, 20] \tag{13}$$

In theory, when the independent variable *x<sup>j</sup>* is equal to 0, the global minimum optimum of the function is *f*(*xj*) = 0. In Figure 5, three different solutions of the function are respectively shown. When the independent variable range is [−20, 20], the function looks smooth and monotonic, as shown in Figure 5a. With the increase of the resolution, many small peaks and valleys are observed in Figure 5b,c. This means that the function has many widespread local minima. ()0 20, 20

( ) 10 cos(2 ) 10 20,20

0

2 0

20, 20 5, 5 0.5, 0.5 **Figure 5.** Rastrigin function: (**a**) *x<sup>i</sup>* ∈ [−20, 20]; (**b**) *x<sup>i</sup>* ∈ [−5, 5]; and (**c**) *x<sup>i</sup>* ∈ [−0.5, 0.5].

To prove the feasibility and accuracy of the modified PSO-MVFSA, two other PSO-based optimization algorithms, as two comparison algorithms, are adopted to search for the optimal minimum of the Rastrigin function. As the function is two-dimensional, each particle consists of two variables for three kinds of PSO algorithms. The initial particle number, the total iteration number of PSO, and the total iteration number of MVFSA are set to 30, 80, and 40, respectively. Some other coefficients are listed in Table 1, such as the positive value α, the Markov chain length *J*, and so on. By means of these searching optimization methods, the optimal results are obtained and shown in Table 2. It can be seen that the optimal function value obtained by PG-PSO is smaller than that of the standard PSO. However, the computing time is longer than that of the standard PSO. For the modified PSO-MVFSA, not only its time consumption, but also its ability for minimum searching, are optimal in comparison with the other two algorithms.


**Table 1.** Coefficients of three algorithms based on particle swarm optimization (PSO).

**Table 2.** Results of three algorithms based on PSO.


#### *3.3. The Improved Airfoil Aerodynamic Optimization*

#### 3.3.1. Fitness Function

For improved cascade aerodynamic optimization, nine control variables are used to parameterize the cascade and grouped as a particle. In all flights, the selection of the particle with the best position through the optimum target is very important. The lift-drag ratio of the airfoil is related to its capability regarding the power output and aerodynamic loss [16]. Therefore, the optimization proposition is set up as shown in Equation (14).

$$\begin{cases} f(\mathbf{S}\_1, \mathbf{S}\_2 \cdots \mathbf{S}\_{10}) = \sum\_{i=1}^n (\mathbf{C}\_L / \mathbf{C}\_D)\_i / n \\ \text{s.t.} \\ \chi\_1 > \arctan \{ B\_y / B\_x \} \\ \alpha\_1 > \arctan \{ T\_y / T\_x \} \\ 0.4 \le D \le 0.6 \end{cases} \quad \begin{aligned} \chi\_1 &\approx \arctan \{ B\_y / (1 - B\_x) \} \\ \alpha\_1 &\approx \arctan \{ T\_y / (1 - T\_x) \} \\ 0.4 \le D \le 0.6 \end{aligned} \tag{14}$$

where *S<sup>i</sup>* is the control variable, *<sup>C</sup>L*/*C<sup>D</sup>* is the lift-drag ratio, *Bx*, *B<sup>y</sup>* is the coordinate of the maximum camber point, *Tx*, *T<sup>y</sup>* is the coordinate of the maximum half-thickness point, α1, α<sup>2</sup> represent the LE and TE thickness gradient angles, χ1, χ<sup>2</sup> represent the LE and TE angles, and *n* is the number of incidence angles. In this paper, the average lift-drag ratio of the airfoil is evaluated under all airflow incidence angles. Additionally, the fitness function can be obtained as shown in Equation (15).

$$\begin{aligned} -F &= c\_1 \frac{f(\mathbf{S}\_1 \mathbf{S}\_2 \cdots \mathbf{S}\_{10}) - \Sigma(\mathbf{C}\_1 / \mathbf{C}\_{\mathrm{D}})\_{\mathrm{ref}}}{\Sigma(\mathbf{C}\_L / \mathbf{C}\_{\mathrm{D}})\_{\mathrm{ref}}} + c\_2 \frac{a\_1 - \arctan\left(\mathcal{B}\_y / \mathcal{B}\_x\right)}{\arctan\left(\mathcal{B}\_y / \mathcal{B}\_x\right)} - c\_3 \frac{a\_2 + \arctan\left(\mathcal{B}\_y / (1 - \mathcal{B}\_x)\right)}{\arctan\left(\mathcal{B}\_y / (1 - \mathcal{B}\_x)\right)} \\ &+ c\_4 \frac{\mathcal{X}\_1 - \arctan\left(\mathcal{T}\_y / \mathcal{T}\_x\right)}{\arctan\left(\mathcal{T}\_y / \mathcal{T}\_x\right)} + c\_5 \frac{-\mathcal{X}\_2 - \arctan\left(\mathcal{T}\_y / (1 - \mathcal{T}\_x)\right)}{\arctan\left(\mathcal{T}\_y / (1 - \mathcal{T}\_x)\right)} + c\_6 (D\_i - D) \end{aligned} \tag{15}$$

where *c<sup>i</sup>* is the weighting factor and subscript *re f* presents the original airfoil.

#### 3.3.2. Aerodynamic Optimization Process

During the evaluation of airfoil aerodynamic performance two-dimensional computational fluid dynamics (2D CFD) code which utilized the standard *k* − ω turbulence model [37,38] to solve the Reynolds Average Navier-Stokes (RANS ) equations was used, and it was only used to simulate the low Reynolds number airfoils [23,39]. Due to the fact that the lift-drag ratio of one airfoil could be obtained quickly by this code, this 2D CFD code was integrated into MATLAB as a fast solver to evaluate the fitness value of each particle. The whole aerodynamic optimization flowchart is shown in Figure 6. Moreover, the detailed process is introduced in the following steps:



 


**Figure 6.** Flow chart of the whole aerodynamic optimization process.

#### **4. Applications in Cascade Optimization**

#### *4.1. Optimization of a Cascade with an NACA4412 Profile*

NACA4412 is adopted as the original basic profile of the cascade and evolved into an advanced airfoil by this improved method. The initial particle number, the total iteration number of the modified PSO, and the total iteration number of MVFSA are set to 30, 100, and 50, respectively. Each cascade is described by the improved aerodynamic parameterization method. The average lift-drag ratio of the airfoil of the cascade is calculated by the flow solver. Moreover, the constant airflow boundary condition is set so that the Reynolds number is as high as 1 × 10 5 , the Mach number is 0.1, and the solidity of the cascade is 1.5. Moreover, some other coefficients of the modified PSO-MVFSA are shown in Table 3.


**Table 3.** Coefficients of modified PSO-modified very fast simulated annealing (MVFSA).

 In actual airflow conditions, the incidence angle of the cascade is not constant and limited to a certain range. Therefore, the aim of multi-point optimization is not to obtain the best aerodynamic performance under one constant incidence angle, but to reach the best whole aerodynamic performance in the whole working range, referred to in Equation (12). In this paper, the incidence angle is uniformly distributed from 0 ◦ to 15 ◦ by a step of 1.0 ◦ . The aerodynamic performance of each cascade parameterized by one incidence angle and eight geometric variables is calculated by the CFD simulation. Based on Equation (13), the fitness value corresponded to each cascade is figured out. Then, the global minimum value is found by the user of the improved PSO-MVFSA method. NACA4412 airfoil and the optimal airfoil are shown in Figure 7, in which the blue curve is the optimal airfoil and the red curve denotes the original airfoil. From this figure, it can be observed that the suction side and pressure side of the optimized airfoil are changed in terms of geometry in comparison with that of NACA4412.

**Figure 7.** Geometry comparison: NACA4412 airfoil and optimal airfoil.

To evaluate the aerodynamic performance of the cascade, the pressure coefficient *C<sup>p</sup>* is adopted, which is defined by Equation (16).

$$\mathcal{C}\_{p} = -\left(\frac{p - p\_{\infty}}{\frac{1}{2}\rho u\_{\infty}^{2}}\right) = \frac{p\_{\infty} - p}{p\_{0} - p\_{\infty}},\tag{16}$$

 where *p*, *p*0, *p*<sup>∞</sup> are the current pressure, the stagnation pressure, and the inlet airflow pressure, respectively; ρ is the density; and *u*<sup>∞</sup> is the inlet airflow velocity.

 *p* The increasing curvature of the suction side close to LE can lead to the intensifying of the velocity increasing and the pressure decreasing. Then, the pressure coefficient of the suction side of the optimized *C<sup>p</sup>* near LE is larger than that of the original cascade. However, it is due to the geometric change of the pressure side close to LE that the decreasing of the velocity and the increasing of the pressure are alleviated. These phenomena can be observed in Figure 8. In the figure, the detailed pressure distributions on the suction and pressure side are shown, which respectively correspond to four incidence angle conditions of *i* = 0 ◦ , *i* = 3 ◦ , *i* = 6 ◦ , and *i* = 12 ◦ . It is also clear that the optimized airfoil surface pressure distribution presented by the blue curve is better than that of the original airfoil described by the red curve under each condition. Therefore, the pressure differences between the pressure side and suction side of the optimized airfoil are larger than those of the original airfoil. The cascade performances under the incidence angles ranging from 0 ◦ to 15 ◦ are shown in Figure 9. The average lift-drag ratio of the optimized cascade is increased by 13.38% in comparison with that of the cascade with NACA4412.

0 3 6 12

°

°

**Figure 8.** Surface pressure distribution comparison in multi-point optimization.

 **Figure 9.** Airfoil performance of the average lift-drag ratio (cal., Re = 1.0 × 10 5 , *Ma* = 0.1).

<sup>5</sup> Re 1.0 10 , 0.1 0 12 For further analysis, CFD simulation software CFX was used to calculate the cascade performance. The velocity contours of the original and optimized cascades at incidence angles of *i* = 0 ◦ and *i* = 12 ◦ are shown in Figure 10. It is clear that the velocity difference between the suction side and pressure side of the optimized cascade is bigger than that of NACA4412. Therefore, the diffusion factor *D* of the optimized cascade can be increased. It can also be observed that the airflow separation point on the suction side of the optimized cascade moves backward along the airfoil in comparison with that of NACA4412 and the area of the wake becomes smaller than that of NACA4412. Therefore, it can be

considered that the aerodynamic load of the cascade is increased and the aerodynamic loss is controlled with the help of the improved aerodynamic optimization algorithm.

**Figure 10.** Velocity contours.

#### *4.2. Blade D500 Optimization*

Blade D500 is used in an axial-flow fan for an evaporator system. It is a kind of low Reynolds number 3D blade. In this work, in order to verify the feasibility of the improved method applied to the practical blade design, the cascade at the 50% radius of Blades D500 was selected and optimized. Two performance parameters, including the pressure coefficient *C<sup>p</sup>* defined by Equation (16) and the efficiency η defined by Equation (17), were used to evaluate the aerodynamic performance of Blade D500. In order to obtain the aerodynamic performance parameters, CFX was used to calculate the airflow field of Blade D500.

$$
\eta = \frac{\mathbb{Q}\_{\mathcal{V}} \cdot P}{N} \,\tag{17}
$$

where *Q<sup>v</sup>* is the volume flow rate, *P* is the static pressure increase, and *N* the aerodynamic power.

#### 4.2.1. Validation of the CFD Simulation Based on Experiments

The experiments of the axial fan with Blade D500 shown in Figure 11 were conducted in the Key Lab for Power Machinery and Engineering of SJTU. In order to match the practical situation, the rotor was mounted on the guard grill, which was connected with a short bell month at the inlet of the tube. The diameter of the tube was increased to 600 mm to avoid destroying the flow field at the rotor outlet and simulate the real situation in the evaporator system as much as possible. The pressure probes were mounted at points A and B, by which the flow rate and the pressure increase could be calculated. Under eight volume flow rates, the pressure increases were calculated.

(**c**) A sketch of the experiment

**Figure 11.** The experiment of Blade D500.

Under the same conditions as those used in the experiments, CFX was adopted to simulate the aerodynamic performance of the axial fan. The grids consisted of the structured hexahedral meshes shown in Figure 12 that were generated by TurboGrid. The total number of nodes was 1,260,626 for the full channel after studying the grid independence, and the minimum value of the mesh quality was 0.3. Due to the excellent ability to simulate the flow with a fiercely adverse pressure gradient of the standard *k* − ω model, it was selected as the turbulent model to simulate the flow field. The mass flow and static pressure were set at the inlet and outlet, respectively.

**Figure 12.** Grid of the axial fan.

After a series of experiments and CFD simulations of the axial fan with Blade D500, the pressure coefficients were obtained and are shown in Figure 13. Based on Figure 13, it could be found that the pressure coefficients *Csp* calculated by CFD are slightly higher than those calculated by experiments. However, the trends of the two curves are very similar. In Figure 14, the efficiency is compared for values obtained by calculations and experiments. It can be observed that the relative error between the

two series of average efficiencies is under 8.35%. Considering the abovementioned results, the CFD simulation can be used to feasibly estimate the aerodynamic performance of an axial fan.

**Figure 13.** Static pressure coefficients calculated by CFD and the experiment.

**Figure 14.** Average efficiencies calculated by CFD and the experiment.

#### 4.2.2. Optimization of Blade D500

 .1 .6 .1 .6 For aerodynamic optimization of the cascade at the 50% radius of Blade D500, the airflow condition was set as Re = 1 × 10 5 , *Ma* = 0.1, *i* = 2 ◦ ∼ 10 ◦ , *D* = 0.4 ∼ 0.6. After a series of optimization iterations, the airfoil of the optimal cascade was obtained and is shown in Figure 15, as well as the airfoil of the original cascade. It can be observed that the airfoil of the optimized cascade is different from that of the original cascade. These geometric changes of the optimal cascade can lead to an increase of the curvature of the suction side, while the pressure side is changed a little. Additionally, aerodynamic performance comparisons are shown in Figure 16. It could be found that the pressure differences of the optimized cascade are larger than those of the original airfoil along the chord line within the range of the incidence angle. After a series of simulation calculations, it could be found that the average lift-drag ratio of the optimized cascade is increased by 15.21% in comparison with that of the original cascade. Therefore, it can be considered that the aerodynamic performance of the optimized cascade is better than that of the original one.

<sup>5</sup> Re 1 10 0.1 2 1 0 0.4 0.6

**Figure 15.** Geometric comparisons of the optimized and original cascades.

**Figure 16.** Surface pressure distributions.

To evaluate the feasibility of the improved aerodynamic optimization method, the optimized 3D Blade D500 and the original 3D Blade D500 were simulated by CFD under the same airflow condition. Moreover, the simulation results of two blades were compared. The streamlines of the suction surfaces of the two blades are shown in Figure 17. From the figure, it is clear that the streamline between the middle section and the hub of the original blade is not very good, and a large turbulence loss must have occurred in this region. In contrast, the streamline in the same region of the optimized blade becomes very smooth. The pressure coefficients *Csp* of the two fans are shown in Figure 18. In this figure, it can be seen that the pressure coefficient of the optimized fan is better than that of the original fan, and the average pressure coefficient is increased by 6.12%. Meanwhile, from Figure 19, it can be observed that the efficiency of the optimized fan is higher than that of the original one, and averagely increased by 11.15%. Therefore, from the above analysis, it can be considered that the aerodynamic performance of optimized Blade D500 is improved and the improved aerodynamic optimization method is feasible.

(**a**) Suction surface of the original cascade (**b**) Suction surface of the optimized cascade

**Figure 17.** Surface streamline comparisons of Blade D500.

**Figure 19.** Efficiency comparisons of D500.

#### **5. Conclusions**

In this paper, an improved optimization method especially applied to a cascade with a low Reynolds number is proposed. In this method, the incidence angle and eight geometric parameters are used as the control variables altogether. Meanwhile, the modified PSO-MVFSA algorithm is adopted to optimize the objective cascades, and two cascade cases, such as NACA4412 and Blade D500, are selected to testify the method. Some valuable results are as follows:


— **Author Contributions:** S.Z. and B.Y., investigation, methodology, software, and writing—original draft preparation. H.X. and M.S., validation. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

Ⅱ **Acknowledgments:** This research was supported by the National Science and Technology Major Project of China (2017-II-0006-0019).

**Conflicts of Interest:** The authors declare that they have no conflict of interest with regards to this work. They declare that they do not have any commercial or associative interests that represent a conflict of interest in connection with the work submitted.

### **Abbreviations**


#### **Nomenclature**


#### **References**


© 2020 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/).

*Article*
