**Experimental and Comparative RANS**/**URANS Investigations on the E**ff**ect of Radius of Volute Tongue on the Aerodynamics and Aeroacoustics of a Sirocco Fan**

#### **Xiaocheng Rui <sup>1</sup> , Limin Lin <sup>2</sup> , Junkui Wang <sup>3</sup> , Xinxue Ye <sup>2</sup> , Haijiang He <sup>2</sup> , Wei Zhang 1,\* and Zuchao Zhu <sup>1</sup>**


Received: 19 October 2020; Accepted: 9 November 2020; Published: 11 November 2020

**Abstract:** The geometry of volute tongue is crucial in the design of Sirocco fans. The size of the volute tongue determines its relative position and distance from the impeller which affects the local flow characteristics and thus the aerodynamic and aeroacoustic performances of the fan. In this work, we performed experimental and numerical investigations on the effect of volute tongue radius on the aerodynamic and aeroacoustic characteristics of a Sirocco fan. The internal flow characteristics are analyzed and discussed in terms of the spatial distribution and temporal variation of pressure and streamlines, the pulsating behaviors of pressure both in the impeller and on the volute surface with emphasis in the volute tongue region, the variation of passage flow with the rotation of impeller and the aeroacoustic features of the fan. We conducted numerical simulations using both steady Reynolds-Averaged Navier-Stokes (RANS) and unsteady Reynolds-Averaged Navier-Stokes (URANS) approaches with realizable *k*-ε turbulence model with rotation effect correction and the results are compared against the experimental data to assess the prediction capability and accuracy in qualitative and quantitative manners. Experimental and numerical results show that as the volute tongue radius increases, the static pressure rises as well as the far-field noise of the fan and pronounced fluctuation of flow is observed within the whole impeller and volute; the reversed flow in the passage of the impeller is reduced and the high-pressure region is found to be moving towards the outlet of the volute. The decreasing radius also enlarges the size of the adverse pressure gradient (APG) region on the volute tongue which contributes to the formation of recirculating flow. The comparative RANS and URANS simulations reveal that both approaches produce generally consistent results regarding the time-averaged flow although the URANS data are much closer to those of the experimental ones. However, the fluctuating flow which is not capable to be modeled by RANS still dominates for the present configuration and thus URANS is necessary for the accurate prediction of the flow details.

**Keywords:** Sirocco fan; URANS; volute tongue radius; internal flow; noise

#### **1. Introduction**

Sirocco fans, also named as multi-blade centrifugal fans, are widely used in industrial and civil applications due to its small size, large flow coefficient and high pressure coefficient, while the

disadvantages of low static pressure and low efficiency have always restricted the development of the Sirocco fans. The analysis and research of the internal flow of the Sirocco fans by means of reasonable and applicable approaches play a meaningful guiding role in further improving and optimizing the performance of the fans.

The volute is one of the main components of the Sirocco fans which plays an important role in guiding the flow outside of the impeller and converting kinetic energy into pressure energy and has received wide extensive attention from many researchers. Rafael et al. applied a numerical model to compute the pressure fluctuation on the volute wall of the fan and obtained an important peak of pressure fluctuation near the volute tongue [1]. By controlling the volute tongue which determines the recirculation flow, Dilin et al. obtained the velocity and pressure distribution under a wide range of azimuth angles [2]. Because the design of the recirculation channel and the volute tongue have a greater impact on the performance and flow stability of the centrifugal fans, Okauchi et al. used four alternative volute designs and accurately measured the flow near the volute outlet and volute tongue and studied the effect of the volute structure on the performance of the fan. The results show that the optimized volute improves the performance of the fan but the flow is unstable when the flow rate increases [3]. Whitfield and Johnson took the compressor as the research object and studied the influence of the volute design on its performance. It was found that while improving the performance of the whole machine, the design method of reducing the flow area under the condition of large flow rate does not have a better effect [4]. Xiao et al. used the energy gradient theory to study the flow instability of the centrifugal fan and finally showed that the impeller outlet and the volute tongue area are prone to produce flow instability and the flow near the hub was more stable than that near the shroud [5]. Pan et al. extended the initial conceptual design of the centrifugal fan and compressor volute and the detailed internal flow data under design and off-design conditions obtained by numerical simulations showed that without changing the ratio of the volute outlet area, the flow area of the volute adjacent to the volute region is increased and the flow distribution is improved [6]. Lun et al. studied the effect of inclined volute on the internal flow and performance of centrifugal fans and found an optimal incline angle to obtain better performance improvement [7]. Sandra et al. conducted an experimental study on the two radial outlet positions of a certain forward curved blade centrifugal fan by hot wire techniques and found a strong flow asymmetry, especially in the circumferential position near the volute tongue and other ones and there is a big difference between the locations [8]. Sina and Javad studied the effect of the volute spread angle on the efficiency and internal flow pattern of Sirocco fans by experimental and numerical methods. The results showed that in a finite space, a fan with a spread angle of 5◦ has the highest head and efficiency [9]. Zhou and Li studied the characteristics of the centrifugal fan volute based on the method of dynamic torque correction to design a better volute profile. It was found that the radial velocity component at the outlet of the volute increases sharply while the noise is reduced without affecting the performance of the fan [10]. Sunil et al. focused on the effect of volute tongue gap on the performance of backward curved blade centrifugal blowers. The comparison of the four sets of data showed that the volute tongue gap has a significant impact on its performance and the total pressure and efficiency of the blower increased with the decrease of the clearance [11].

The impeller is the main component of the Sirocco fans and its geometry, size, number of blades and so forth have a crucial influence on the performance of the fan. Taking a widely used centrifugal fan as the research object, Zhang et al. analyzed the influence of the optimized combination of impeller geometrical parameters, such as the number of blades and the impeller outlet width, on the performance and noise of the fan. Under the design conditions, the combination of parameters after optimization makes the total pressure significantly increasing by 6.91% [12]. Lee et al. proposed a method of redesigning the centrifugal impeller and its intake passage. The resulting effect not only satisfies the pressure requirements but also reduces the power of the fan by 8.8% compared with the baseline model [13]. Wang et al. proposed an inverse design method for calculating the meridional surface of the centrifugal impeller. The final computation results show that the performance of the fan inlet is improved and the velocity distribution is more uniform [14]. Kim et al. conducted a

parameter optimization study on the impeller geometry of the forward curved blade centrifugal fan. The height of the annular plate separating the impeller and the angle between the upper and lower impellers were chosen as geometric parameters. The influence of these parameters on the efficiency of the fan was analyzed. The results show that reasonable parameters are conducive to improving the aerodynamic performance of the fan [15]. Ni et al. changed the straight blades of the baseline Sirocco fan impeller into inclined blades with a certain angle, obtained an optimal value of the inclination angle and improved the aerodynamic performance of the fan [16]. Kim and Seo successfully applied the response surface optimization method of three-dimensional N-S analysis to the aerodynamic design of forward curved blade centrifugal fans and obtained the maximum static pressure efficiency [17].

The influential factors, such as the interaction or positional relationship between the impeller and the volute or other components, have a significant effect on the internal flow of the fans. Tarek and Seung conducted a numerical simulation analysis of the interactions among the impeller, diffuser and volute. It was found that the vortex flow area in the volute is gradually attenuated from the part near the volute tongue to the fan outlet [18]. Koen and Rene proposed a numerical method for predicting the interaction between impeller and volute in a single-stage centrifugal compressor and analyzed the influence of the interaction between impeller and volute on the circumferential distortion of static pressure [19]. Karanth and Sharma considered the complex flow between the impeller outlet and the diffuser, studied and predicted the flow characteristics caused by the radial clearance by numerical simulation analysis and the results showed that there is an optimal radial clearance which made better energy conversion of the impeller [20]. Lee conducted a study and analysis of the gap effect between the inlet duct and the rotating impeller of the centrifugal fan. The results showed that improper gap determination would reduce the overall performance of the fan by 2–5% [21]. Jose et al. studied the interaction between the impeller and the volute and numerically analyzed the secondary flow in the volute through the helicity magnitude. The results showed that there is a strong secondary flow in the radial position near the outlet of the impeller [22]. Daniel et al. considered the influence of the relative position of the impeller and the volute on the performance of the centrifugal pump. Through the research, it was found that the efficiency of circular volute and spiral volute increases by 5% and 3.5%, respectively, as the impeller is in the optimal position compared with the central position [23]. Younsi et al. studied the interaction and unsteady characteristics of the rotating blades relative to the movement of the volute. Analysis of the pressure fluctuations on the volute surface showed that the main source of flow disturbance and instability in the centrifugal fan is the volute region [24].

Bechara et al. devised a stochastic noise generation and radiation (SNGR) model for noise generation and radiation from the turbulent flow field [25]. Billson et al. used a SNGR method to capture the sound field of a high Reynolds number, high Mach number subsonic jet. The directivity of the far-field sound is well predicted but the amplitude and spectral information differ from the measurements [26]. Ewert studied the application of a low-cost computational aeroacoustics (CAA) approach to a slat noise problem. Slat noise simulations are carried out for a high-lift airfoil and the Mach number scaling law of the broadband slat noise component is evaluated based on three different free-stream velocities (M = 0.088, 0.118, 0.165) [27]. In the following work, the simulated 1/3-octave spectrum exhibits the main features of an empirical slat-noise mode [28]. Bailly et al. introduced an application of direct calculation of aerodynamic noise in subsonic and supersonic jet noise, cavity noise and self-excited internal flows [29]. Sandberg conducted direct numerical simulations (DNS) to airfoil self-noise. It is found that the contribution of trailing edge noise dominates at low frequencies while for the high frequencies the radiated noise is mainly due to flow events in the reattachment region on the suction side [30]. Dawi and Akkermans studied spurious noise generated in direct noise computation using finite volume methods. Different sources of spurious noise are examined as well as the mechanism of their generation [31]. In another work, they used direct noise computation and Fflowcs Williams-Hawkings equation to calculate far-field spectral. The peaks of vortex shedding frequency were well predicted, however, the amplitude of the spectra at these frequencies is smaller than the ones in the wind tunnel [32]. Deuse and Sandberg studied the self-noise of an isolated

controlled-diffusion aerofoil using direct noise computation. Two main noise sources were observed, one at the trailing edge and one in the transition/reattachment region. For the cases investigated, the latter is even more important than the trailing edge noise component [33]. Dawi and Akkermans demonstrated the applicability of a finite volume method for the direct noise computation of road vehicles, the results of the DNC showed good agreement with wind tunnel measurements within a certain frequency range. However, at higher frequencies, simulated acceleration spectra showed lower spectral levels than the measured ones [34]. Dierke et al. investigated the effect of installation on propeller sound. It is found that the installation effect due to the presence of the high-lift wing increases the sound pressure levels by 5–10 dB in most directions [35]. The research by Dobrzynski et al. on airframe noise prediction and reduction methods pointed out that the noise level of gears installed under the wings is reduced by about 4 decibels relative to gears installed on the fuselage. The aerodynamic optimum design of current slotted slats causes maximum gap-flow velocities and thus high slat noise levels [36]. Ewert et al. analyzed the changes of acoustic level due to a change of slat position and found that by changing the slat gap, the acoustic level of some frequency segments also changed correspondingly [37].

It can be concluded from the above literature review that the geometry of the volute, especially the volute tongue radius, the distance to the impeller outlet and the axial change of the volute surface determine the flow in the vicinity of the impeller outlet and volute tongue and affect the interaction between the fluid within the impeller and the volute, with a vital influence on the performance of the Sirocco fan. Considering that the radius of the volute tongue determines the distance between the local volute wall and the impeller, in this work we performed experimental and numerical investigation on the internal flow of a Sirocco fan with various volute tongue radii to identify its effect on the aerodynamic and aeroacoustic characteristics of the fan. Only the radius of the volute tongue is varied but the distance to the outlet of the impeller keeps unchanged and the internal flow is analyzed in terms of the spatial distribution and temporal variation of pressure and streamlines, the pulsating behaviors of pressure both in the impeller and on the volute surface with emphasis in the volute tongue region, the variation of passage flow with the rotation of impeller and the aeroacoustic features of the fan. Another objective of the present work is to explore the capability, reliability and accuracy of RANS and URANS approaches in modeling flow in similar fan configurations both in qualitative and quantitative manners and provide meaning conclusions in further simulations.

The paper is organized as follows—Section 2 is a general description of the fan configuration, the employed numerical methods, the experiment details and the validation. Section 3 presents and discusses the numerical and experimental results for the effect of radius of volute tongue. Some conclusions are drawn in Section 4.

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

#### *2.1. Fan Configuration*

The fan in this investigation is a double-suction, forward-curved multi-blades Sirocco fan which is shown in Figure 1. The blades are staggered on both sides of the impeller's central disk with 41 blades on each side. The air enters the volute and impeller from the inlet domain on both sides of the fan and moves out from the outlet domain. Table 1 lists the specifications of the baseline model. The radius of volute tongue of the baseline model is 11 mm and the location is denoted in the figure. In this work, we investigate the effect of volute tongue radius at four values, that is, *r* = 4 mm, 9 mm, 11 mm and 13 mm, on the aerodynamic and aeroacoustic characteristics of the fan at the designed flow rate 365 m<sup>3</sup> /h.

**Figure 1.** Baseline model. (**a**) fan model; (**b**) impeller and volute; (**c**) definition of radius of volute tongue.


**Table 1.** Specifications for the baseline Sirocco fan.

#### *2.2. Numerical Simulation Method*

The internal flow of the Sirocco fan is governed by the conservation of mass and momentum equations:

$$\frac{\partial \mu\_i}{\partial \mathbf{x}\_i} = \mathbf{0} \tag{1}$$

$$\frac{\partial u\_i}{\partial t} + u\_j \frac{\partial (u\_i)}{\partial \mathbf{x}\_j} = f\_i - \frac{1}{\rho} \frac{\partial \mathbf{P^\*}}{\partial \mathbf{x}\_i} + \nu \frac{\partial^2 u\_i}{\partial \mathbf{x}\_j \partial \mathbf{x}\_j} \tag{2}$$

డ௨ ݑ + డ௧ డሺ௨ ሻ డ௫ೕ = ݂ − ଵ ఘ డ<sup>∗</sup> డ௫ డ ߥ + <sup>మ</sup>௨ డ௫ೕడ௫ೕ ݔ ߩ where ρ represents the density of the air; *x<sup>i</sup>* is the direction component of the Cartesian coordinate system; *u<sup>i</sup>* is the velocity component; *P* ∗ is the pressure considering the conversion of turbulent kinetic energy and centrifugal force; *f<sup>i</sup>* is the component of the volume force. ν is the kinematic viscosity coefficient:

$$\nu = \frac{1}{\rho}(\mu + \mu\_t) = \frac{\mu}{\rho} + \mathcal{C}\_{\mu}\frac{k^2}{\varepsilon},\tag{3}$$

where µ is the molecular viscosity coefficient and µ*<sup>t</sup>* is the vortex viscosity coefficient. In the numerical simulation, we use the realizable *k*-ε turbulence model with rotation effect correction in which *C*<sup>µ</sup> is calculated by the following formula:

$$\mathcal{C}\_{\mu} = \frac{1}{A\_0 + A\_S \frac{kL^\*}{\varepsilon}} \,' \tag{4}$$

in which *A*<sup>0</sup> = 4.04, *A<sup>S</sup>* = √ 6 cosφ, φ = <sup>1</sup> 3 cos−<sup>1</sup> √ 6*W* , *W* = *SijSjkSki* e*S* , e*S* = p *SijSij*, *Sij* = <sup>1</sup> 2 ∂*u<sup>j</sup>* ∂*x<sup>i</sup>* + ∂*u<sup>i</sup>* ∂*x<sup>j</sup>* , *U*<sup>∗</sup> ≡ q *<sup>S</sup>ijSij* + <sup>Ω</sup>e*ij* <sup>Ω</sup>e*ij*, <sup>Ω</sup>e*ij* = <sup>Ω</sup>*ij* <sup>−</sup> 2ε*ijk*ω*<sup>k</sup>* , <sup>Ω</sup>*ij* = <sup>Ω</sup>*ij* <sup>−</sup> <sup>ε</sup>*ijk*ω*<sup>k</sup>* . Ω*ij* is the laminar curl with angular velocity.

The transport equations for the turbulent kinetic energy *k* and turbulent dissipation rate ε in realizable *k*-ε turbulence models are as follows:

$$\frac{\partial}{\partial t}(k) + \frac{\partial}{\partial \mathbf{x}\_i}(k\boldsymbol{u}\_j) = \frac{1}{\rho} \frac{\partial}{\partial \mathbf{x}\_i} \left[ \left(\mu + \frac{\mu\_l}{\sigma\_k}\right) \frac{\partial k}{\partial \mathbf{x}\_j} \right] + \mathbf{G}\_k + \mathbf{G}\_b - \rho \boldsymbol{\varepsilon} - \mathbf{Y}\_M + \mathbf{S}\_k \tag{5}$$

$$\frac{\partial}{\partial t}(\varepsilon) + \frac{\partial}{\partial \mathbf{x}\_{j}} (\varepsilon u\_{j}) = \frac{1}{\rho} \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \left( \mu + \frac{\mu\_{l}}{\sigma\_{k}} \right) \frac{\partial \varepsilon}{\partial \mathbf{x}\_{j}} \right] + \mathbf{C}\_{1} \mathbf{S}\_{\varepsilon} - \mathbf{C}\_{2} \frac{\varepsilon^{2}}{k + \sqrt{\nu \varepsilon}} + \mathbf{C}\_{1\varepsilon} \frac{\varepsilon}{k} \mathbf{C}\_{3\varepsilon} \mathbf{G}\_{\theta} + \mathbf{S}\_{\varepsilon}. \tag{6}$$

In the equation, *G<sup>k</sup>* and *G<sup>b</sup>* represent the turbulent kinetic energy terms generated by laminar velocity gradient and buoyancy, respectively. *Y<sup>M</sup>* represents the contribution to the dissipation rate of the expansion of turbulent pulsation to the global flow in compressible flows. *C*1<sup>ε</sup> and *C*<sup>2</sup> are constants. σ*<sup>k</sup>* and σ<sup>ε</sup> are the Prandtl numbers. *S<sup>k</sup>* and *S*<sup>ε</sup> are user-defined turbulent kinetic energy and turbulent dissipation source terms. For the incompressible air as the working medium, *G<sup>b</sup>* = 0 and *Y<sup>M</sup>* = 0.

In this work, we performed both RANS and URANS simulations using ANSYS-Fluent 17.0 (Version, ANSYS, Inc., Canonsburg, PA, USA, 2016) for the same models and flow parameters. In the simulations, second-order upwind scheme is used for the spatial discretization of all derivative terms; the SIMPLE method is used for the coupling of velocity and pressure and the moving reference frame (MRF) method is used to treat the rotation of the impeller. The boundary conditions are the same for both simulation approaches. Constant mass flow rate is prescribed at the two inlet sections and zero gauge pressure at the outlet section. No-slip boundary condition is set for velocity on all solid walls. We first perform the RANS simulation until the converged result with relative residual less than 10−<sup>3</sup> is obtained for the equations of mass, momentum, turbulent kinetic energy and turbulent dissipation rate. The result is then used as initial condition to start the URANS simulation and the simulation continues as the impeller rotates for ten revolutions before the data are collected. The physical time step size is fixed as the impeller rotates for one degree, thus 360 steps are required for one revolution of the impeller.

#### *2.3. Grid Independence Study*

The computational domains of the fan include two inflow domains on both sides, an outflow domain, a volute domain, an impeller domain and the intermediate flow channel domain of the impeller. The impeller and intermediate flow channel of the impeller constitute the rotational domain, while the other domains are set stationary. The impeller, inflow and outflow domains are discretized by block-structured grid, while the volute domain is discretized by unstructured grids considering its complex geometry and subsequent geometrical modifications. Figure 2 shows the grid distributions adjacent to the impeller and volute and the grids number is shown in Table 2.

**Figure 2.** (**a**): The grids around the blade, (**b**): The grids around the volute tongue.


**Table 2.** Number of grids of the computational domain.

In this paper, the static pressure rise, computed as the difference between the static pressure at the outflow and inflow sections of the fan, is used as the indicator for the grid independence study under design flow rate of 365 m<sup>3</sup> /h. The grid number is changed mainly for the impeller domain. Three sets of impeller grids with the number of about 4.6 million, 5.7 million and 6.7 million are used, respectively. The variation of static pressure rise with grid number is shown in Figure 3. As the number of grids increases from 4.6 million to 5.7 million, the variation of static pressure is 1.07 Pa; as the grid number continues to change to 6.7 million, the variation of static pressure rise is only 0.08 Pa, which is only 0.1% higher than the static pressure at 5.7 million. Therefore, we consider that the influence of a higher grid number on the computational results can be negligible, so the impeller domain with a grid number about 5.7 million is chosen for subsequent work.

**Figure 3.** Variation of static pressure rise with grid number of impeller domain.

#### *2.4. Experiment Details*

The various volute models of the Sirocco fan are made by 3D printing and are experimentally tested in the aerodynamic-aeroacoustic laboratory of Zhejiang Yilida Ventilator Co., Ltd. (Taizhou, China). The experimental device is designed according to National Standard of China GB/T1236-2017 and conforms to the International Standard ISO5801-2007. The semi-anechoic noise test is designed according to ISO3745-2003. The noise monitoring can be carried out synchronously when the aerodynamic performance of the fan is tested. The volute models and the impeller are comparably shown in Figure 4. The schematic diagram and photo of the experimental facility are shown in Figure 5. The size of the east and west semi-anechoic room is 5200 mm × 4800 mm × 7000 mm and 5200 mm × 5600 mm × 7000 mm, respectively and the air compartment size is 3000 mm × 3000 mm × 9000

mm. The semi-anechoic room and the air compartment are connected by a circular pipe and the overall loop is a closed-chamber circulation structure. The measurement system mainly includes fan, motor, photoelectric infrared speed sensor, microphone and so forth, as shown in Figure 6. The test facility could measure and record the pressure and temperature of the inlet and outlet flows and the aerodynamic noise in a synchronous way.

**Figure 4.** Impeller and volute with tongue radius *r* = 4 mm, 9 mm, 11 mm and 13 mm.

**Figure 5.** (**a**) The schematic diagram, (**b**) photo of the facility.

**Figure 6.** Test equipment: (**a**) fan and motor; (**b**) photoelectric infrared speed sensor; (**c**) microphone.

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

#### *3.1. Aerodynamic Performances*

The static pressure rise and static pressure efficiency are the key parameters reflecting the aerodynamic performances of the Sirocco fan and is experimentally measured as well as simulated in this work. Figure 7 gives the time history and time-averaged value of the static pressure rise of the fan obtained by URANS and compared with those obtained by RANS simulation and experiment. Tables 3 and 4 list the static pressure rise and efficiency predicted by experiment, RANS and URANS in time-averaged sense.

**Figure 7.** Time history of static pressure and its averaged value computed by unsteady Reynolds-Averaged Navier-Stokes (URANS) in comparison with those of Reynolds-Averaged Navier-Stokes (RANS) and experiment.


**Table 3.** Experimental and numerical results of static pressure.



It is seen from the figure and tables that for models with *r* = 9 mm, 11 mm and 13 mm, the relative difference between the results of RANS and experiment is 6.3% at maximum which is of engineering significance. However, the URANS produce much better consistent data with respect to the experimental one; the maximum relative difference between the time-averaged static pressure rise and experimental data is only 0.7%. This reflects that the URANS simulation is more consistent with the realistic operation of the fan and is capable to simulate the flow inside the fan more accurately. For the model with *r* = 4 mm, the sharper volute tongue strongly destabilize the local flow around the outlet of the volute, thus the pressure rise substantially fluctuates with a large amplitude and the periodicity due to impeller rotation is clearly observed. Both the results produced by RANS and

URANS in time-averaged sense notably deviate from the experimental data with relative difference about 10%. Therefore, we believe that the accuracy of simulation is greatly affected by the pulsating flow and the simulation is less capable in accurately predicting the internal flow of the fan.

Comparing the models with different volute tongue radius, it can be found that the pressure fluctuation in URANS from the 5th to the 10th revolution of the impeller presents the smallest amplitude for the *r* = 13 mm volute, which almost perfectly coincides with the time-average value. For the *r* = 9 mm and 11 mm models, the amplitude slightly increases but is still small compared with the magnitude of pressure rise. The amplitude is the highest for the *r* = 4 mm model. It is also noted that as the radius decreases, the time-averaged pressure monotonically increases from *r* = 13 mm to *r* = 9 mm, while the increasing generally stops for *r* = 4 mm. This indicates that smaller radius will generate stronger pressure pulsation and increases the pressure rise to different degrees. For the *r* = 4 mm model, notable quasi-periodic pressure fluctuation is observed in which one cycle approximately corresponds to 0.05 s as the revolution of the impeller; this is an indication that the sharper volute tongue could generate strong pulsating flow compared with other models.

#### *3.2. General Characteristics of Internal Flow*

The transient flow in the impeller and volute is analyzed for the 10th revolution of the impeller at four moments, that is, T1, T2, T3 and T4, with time interval between neighboring two moments as the impeller rotates for 90◦ . Figures 8–11 show the distribution of static pressure on axial cross-section Z = 4 mm which is close to the central disc of the impeller.

**Figure 8.** Distribution of static pressure on Z = 4 mm cross-section for the *r* = 4 mm model.

**Figure 9.** Distribution of static pressure on Z = 4 mm cross-section for the *r* = 9 mm model.

**Figure 10.** Distribution of static pressure on Z = 4 mm cross-section for the *r* = 11 mm model.

**Figure 11.** Distribution of static pressure on Z = 4 mm cross-section for the *r* = 13 mm model.

It is clearly seen that for models with *r* = 9 mm, 11 mm and 13 mm, the pressure distribution is quite similar at different moments in both the impeller and volute which is consistent with the conclusion obtained in Figure 7. The pressure has a minor fluctuation not only exhibited by its global behavior as the static pressure rise of the fan but also in a local manner. For the *r* = 4 mm model, noticeable temporal variations of pressure field are observed especially in the region between the volute tongue and outlet of the impeller. At time T2, the high-pressure area near the volute tongue is substantially smaller than that at the moment T1, while it increases in size at T3 and continues increasing to a certain extent at T4. Meanwhile, the low-pressure region that covers most of the impeller at time T2 has a wider distribution and a lower pressure magnitude. The above phenomenon indicates that the velocity of passage flow in the impeller is higher at T2 and the flow moves smoothly around the volute tongue. However, the flow is greatly affected by the sharp volute tongue and experiences pronounced changes.

Comparing the pressure distribution of different models, it can be found that for the *r* = 9 mm, 11 mm and 13 mm models, the high-pressure region in the vicinity of the volute tongue slightly expands in area, while it is notably smaller for the *r* = 4 mm model. As the flow velocity is in a state of change, at T2, part of the impeller passage exhibits a certain range of low pressure areas due to the higher flow velocity. This indicates that the volute tongue radius not only affects the flow near the volute tongue but also indirectly affects the flow in other region away from it.

#### *3.3. Reversed Flow in Impeller Near the Volute Tongue*

#### 3.3.1. Reversed Flow Predicted by RANS

It is seen in Figures 8–11 that there is great pressure fluctuation in the passages of the impeller near the volute tongue, thus reversed flow, that is, in the direction from the outlet to inlet of the impeller, may form during part of the revolution cycle. In order to explore the influence of volute tongue radius on reversed flow, in this section we present and analyze the flow in the impeller passages near the volute tongue. Figure 12 gives the schematic of the eight selected impeller passages. The angular cross-section at the 50% spanwise height is chosen as a representation covering the passages around the volute tongue and is expanded as on the X-Z plane in the figure. The distribution of *v*-velocity of flow in the selected passages is plotted in Figure 13. Since the selected blade passages are basically located in the region between 30◦ and 100◦ , the positive and negative *v*-velocity could generally indicate the forward passage flow which moves from the inlet to outlet of the blade passage and the reversed flow moving from the outlet to the inlet of the impeller.

**Figure 12.** (**a**) Selected blade passages near the volute tongue, (**b**) the expanded view of the passages.

It is seen in Figure 13 that there is almost no reversed flow in passage-7 and passage-8 as predicted by the steady-state RANS simulation, while reversed flow of different scales and intensity form in passages 1–6. There is a large region with strong reversed flow close to the suction surface of the impeller and the central disc of the impeller, while forward flow is mainly observed near the pressure surface of the blade. Reversed flow starts to form in passage-6; the reversed flow with small magnitude occupies only a small region for the *r* = 4 mm mode but gradually expands in size and intensified as the radius increases. The same variation pattern of the reversed flow is also observed in passage-5. This variation pattern is attributed to the fact that the volute surface covers almost the whole passage-5 and passage-6 for the *r* = 4 mm model that the passage flow is not affected by the volute tongue, while for the models of *r* = 9–13 mm the volute tongue locates right beside the two passages, thus the flow is influenced by the local flow in the volute and presents reversed flow. For flow in passage 1–4, it is found that as the radius increases, there is more forward flow in the corner region of the pressure surface of the blade and the front disc of the impeller and the reversed flow still exists especially for passage-3 and passage-4. In general, the results of steady-state RANS simulation could reflect the primary characteristics of flow direction in the blade passages near the volute tongue.

**Figure 13.** Distribution of *v*-velocity in the blade passages near the volute tongue.

#### 3.3.2. Reversed Flow Predicted by URANS

In this section, the temporal variation of passage flow in the impeller is analyzed based on the transient results obtained by URANS simulation. We choose three fixed points in passage-1 at the middle between the pressure and suction side of the blade and on the Z = 4 mm, 45 mm and 90 mm axial cross-sections as the monitoring points and record the temporal variation of radial (*r*-) velocity as the passage moves to a series of circumferential positions with the rotation of impeller, as shown in Figure 14a, thus the variation of velocity reflect the temporal variation of the passage flow. The radial velocity is computed by vector manipulation based on *u*-velocity and *v*-velocity and forward and backward flows are identified as the radial velocity is positive and negative, respectively. To focus on the effect of volute tongue on the passage flow in the impeller, more circumferential positions are selected in the region near the volute tongue, while fewer monitoring points are positioned in the rest part of the circumference. Moreover, to facilitate the direct comparison with the steady-state RANS result, we also computed the radial velocity based on the result of RANS simulation at several circumferential positions which are identical with those for the URANS result, as shown in Figure 14b. Since the impeller is axisymmetric about the center of the geometry, the RANS result presents the variation of radial velocity under the Euler framework, while the URANS result is obtained under the Lagrange framework. The difference between the RANS and URANS results reveal the different capabilities of the two approaches.

Figure 15 presents the circumferential distribution of radial velocity for monitoring points on the Z = 4 mm, 45 mm and 90 mm axial cross-sections. For the Z = 4 mm cross-section close to the central disc of the double-suction impeller, the transient results show that there are primarily three regions for all four models categorized by the radial velocity, that is, the forward flow region roughly lower than 30◦ , the reversed flow region in the range between 30◦ and 130◦ and forward flow region after 130◦ . For the *r* = 4 mm model, the reversed flow is observed to occupy only a small region and the magnitude of radial velocity is lower compared with other three models and the flow turns to be forward around

70◦ . As the volute tongue radius increases to *r* = 13 mm, stronger reversed flow generates within this region. Referring Figure 14a, it can be concluded that as the blade passage approaching to the circumferential position Φ1, the magnitude of radial velocity of the reversed flow first increases and then decreases and then exhibits a local maximum near Φ2. The reduction of the volute tongue radius could reduce the range and velocity magnitude of reversed flow. In the range 0–30◦ , the radial velocity gradually decreases until it reaches zero and gradually increases in the range 130–360◦ . Comparing the RANS and URANS results, we can notice that there is great difference between the curves for both forward and reversed flows and for almost the entire circumference of the impeller which may be attributed to the existence of the central disc. For the Z = 45 mm cross-section, the transient result show that reversed flow also exists in quite small regions around the circumferential positions 30◦ and 120◦ , that is, close to Φ1 and Φ2 as shown in Figure 14a. The variation pattern of radial velocity in the circumferential direction is similar with that of the Z = 4 mm cross-section. A comparison of RANS and URANS results demonstrate that significant differences in the range 320–400◦ (40◦ ) and 80–180◦ ; in particular, the flow can be falsely predicted as reversed flow by RANS and the magnitude of velocity also deviate a lot compared with the URANS data. For the Z = 90 mm cross-section close to the front disc of the impeller, the transient result shows that reversed flow appear in the regions 0–40◦ , 110–130◦ and 200–270◦ although the magnitude of radial velocity is minor. The influence of volute tongue radius on the radial velocity is comparably smaller than those on the Z = 4 mm and 45 mm cross-sections. Pronounced difference is still observed between the RANS and URANS results.

**Figure 14.** (**a**) Positions for the transient motion of the fixed point where pressure is recorded in URANS simulation; (**b**) positions where pressure is recorded in RANS simulation.

Since the flow enters the impeller from the front disc, the velocity direction changes from axial to radial as it gradually approaches the central disc. Therefore, near the front disc of the impeller, the radial velocity of the flow is relatively small and the degree of velocity change is small, which can be observed by comparing the distribution of the radial velocity on the three sections. The reversed flow roughly in the region 30–130◦ has the smallest and largest range on the Z = 45 mm and Z = 90 mm cross-sections in RANS and URANS results, respectively, which reflects that the flow near the central disc is less affected by the rotation of impeller. The comparison of RANS and URANS results show that the circumferential variation of velocity magnitude is much larger for the RANS simulation particularly on the Z = 4 mm cross-section where the RANS result is also greatly affected by the volute tongue radius.

Φ Φ

Φ

Φ

**Figure 15.** Variation of radial velocity within the blade passages.

#### *3.4. Near-Wall Flow of the Volute Surface*

For the Sirocco fan investigated in this work, the geometry of the volute is modified only for its tongue, while the rest majority part of the volute surface keeps unchanged. In this section, we analyze the near-wall flow of the volute surface to explore how the modification on the volute tongue affects the flow around the whole circumference of the volute. The near-wall flow is studied on the Z = 4 mm axial cross-section close to the central disc at four moments T1, T2, T3 and T4 same as those shown in Figures 8–11 and the steady-state RANS result is also provided for comparison. The distributions of pressure coefficient and skin friction coefficient will be given on a local coordinate S along the surface of the volute with the origin chosen at the exit of the volute, as shown in Figure 14a.

#### 3.4.1. Distribution of Skin Friction Coefficient on the Volute Surface

Figure 16 shows the distribution of skin friction coefficient (*C<sup>f</sup>* ) around the whole surface of the volute (left subfigure) and the enlarge view for the curves around the volute tongue (right subfigure). For the *r* = 4 mm model, the magnitude of *C<sup>f</sup>* drastically increases in the volute tongue region and the local near-wall flow is almost steady-state from S = 0 to S = Φ1, as reflected by the overlapping curves of various moments; this indicates that although the flow in the impeller is highly fluctuating (as seen in Figure 8), it hardly impose obvious fluctuation for the near-wall flow of the volute. we have observed two peaks of *C<sup>f</sup>* and one valley in the region around the volute tongue and the magnitude of *C<sup>f</sup>* changes dramatically, This is because it is not only close to the outlet of the volute but also has significant changes in wall structure, so the flow state here is the most complex and diverse. The temporal fluctuation of *C<sup>f</sup>* in the rest part of the volute surface is relatively significant. In the region from Φ1 to Φ2, the magnitude of *C<sup>f</sup>* gradually decreases which indicates the decelerated flow. From Φ2 to Φ3, the magnitude of *C<sup>f</sup>* increases first and then decreases. Since Φ2 and Φ3 are located at

the two corners of the bottom of the volute, it suggests that the abrupt change of volute geometry will affect the local near-wall flow. From Φ3 to Φ5, the magnitude of *C<sup>f</sup>* first increases rapidly and then slowly decreases to varying degrees. This variation pattern may be caused by the varied geometry of volute surface from a flat one to a curved one and then gradually varies as an arc and it is also related to the variation of spacing between the impeller and volute surface. In general, in the majority part from Φ3 to Φ5, the variation of *C<sup>f</sup>* is mild as well as the near-wall flow. For the surface at the outlet section of the volute beyond Φ5, the distance between the volute and the impeller is large that the local flow is decelerated, thus the magnitude of *C<sup>f</sup>* also reduces.

For other models of *r* = 9 mm, 11 mm and 13 mm, the curves of *C<sup>f</sup>* shows a certain degree of temporal fluctuations and the variation pattern is basically consistent with that of *r* = 4 mm. Comparing the several models, it can be observed that in the range between S = 0 and S = Φ1, there are always two local peaks. The difference is that with decreasing radius, the slope of *C<sup>f</sup>* curve before and after the valley point is higher and the variation is more intense, which indicates that the decreasing radius generates a smaller low-velocity region where the flow is less affected and thus is beneficial to the aerodynamic performance of the fan. Moreover, the magnitude of second peak reduces with increasing radius from about 5.7 at *r* = 4 mm to 4.7 at *r* = 13 mm; it means small radius is more likely to produce a high-velocity flow at the initial position of the recirculating flow at the tongue. In the region beyond Φ1, the unsteadiness of near-wall flow is observed to be quite different on the volute surface; substantial fluctuation is found beyond Φ1, Φ3, Φ2 and Φ1 for the four models but the fluctuation is generally the strongest for the *r* = 4 mm model. The RANS and URANS simulations have generally consistent predictions on the variation of *C<sup>f</sup>* around the volute surface especially in the region beyond Φ1. There is a significant gap between the results of two approaches in the region from S = 0 to S = Φ1 and the RANS approach fails to capture the two peaks, which demonstrates that although RANS could only depict the general variation of near-wall flow but not capable for the details in the region with great geometrical variation.

#### 3.4.2. Distribution of Pressure Coefficient on the Volute Surface

The near-wall flow of the volute surface is greatly affected by the pressure gradient field. Figure 17 shows the distribution of pressure coefficient (*Cp*) around the whole surface of the volute (left subfigure) and the enlarge view for the curves around the volute tongue (right subfigure). For the *r* = 4 mm model, the distribution of *C<sup>p</sup>* varies generally as a whole in time, that is, the magnitude of *C<sup>p</sup>* at different position synchronously increases or decreases, which indicates that the pressure on the volute surface has obvious unsteady characteristics during the whole cycle of the revolution of impeller. In the region from S = 0 to S = Φ1, there is noticeable fluctuation of *C<sup>p</sup>* compared with the quasi-steady *C<sup>f</sup>* , while the pressure gradient field, as measured by the slope of the curve, does not vary a lot during the whole cycle, thus the velocity of near-wall flow does not show remarkable temporal variation. From the peak value at S = 0, favorable pressure gradient (FPG) with a sharp decrease and then a relatively slow decrease of pressure deteriorates the movement of flow towards the outlet of volute. The APG field from the valley to Φ1 will accelerates the near-wall flow. At the volute tongue, the sharp increase and decrease of *C<sup>p</sup>* clearly show the pressure variation characteristics under the impact of the flow. There is an adverse pressure gradient field (*dp*/*ds* > 0) from Φ1 to Φ2 which is detrimental to the movement of the near-wall flow and may even cause flow separation. From Φ2 to Φ3, the pressure coefficient first decreases and then increases, indicating that the flow velocity may be accelerated and then decelerated in this region. In the region from Φ3 to Φ5, the magnitude of *C<sup>p</sup>* varies relatively smoothly with only minor increase. It is found that the distance between the volute surface and the impeller gradually expands within this region, the geometry of the volute surface abruptly changes from a flat surface to a curved surface and then slowly changes to an arc, which makes the flow smooth in this region. In the region after Φ5, the magnitude of *C<sup>p</sup>* decreases and then increases again until the outlet of the volute because the distance between the impeller and the volute surface increases rapidly which reduces the pressure. For other three models, the variation of *C<sup>p</sup>* is similar to that of the *r* = 4 mm model.

**Figure 16.** Distribution of skin friction coefficient around the volute surface (left column) and the enlarged view around the volute tongue (right column) on the Z = 4 mm cross-section.

(**d**) *r* = 13 mm

**Figure 17.** Distribution of pressure coefficient around the volute surface (left column) and the enlarged view around the volute tongue (right column) on the Z = 4 mm cross-section.

Comparing the *C<sup>p</sup>* distribution of four models, it can be seen that the variation is quite similar for the *r* = 9–13 mm models where the temporal variation is negligible, while the highly fluctuating flow induce noticeable pressure fluctuation for the *r* = 4 mm model which is consistent with the observations in Figure 8. Around the volute tongue, the *C<sup>p</sup>* curve has the steepest increasing and decreasing trend, while the variation is less pronounced for other three models.

The variation of *C<sup>p</sup>* has similar pattern as predicted by RANS and URANS approaches although the magnitude deviates to some degree; the steady-state RANS result could generally reflect the variation but the magnitude is still not accurately predicted even for the *r* = 9–13 mm models where the flow is quasi-steady.

#### *3.5. Staic Pressure Fluctuation and Its Propagation*

From the discussions above, we can conclude that the fluctuation of flow is dependent on the volute tongue radius and is different in the impeller or in the volute. To further analyze the static pressure fluctuation, a series of static pressure monitoring points are set on the Z = 4 mm axial cross-section near the central disc where the time history of pressure are recorded during the 8th to 10th revolution in the URANS simulation, as shown in Figure 18. The monitoring points are located in the region where the temporal variation of local flow is significant and can be compared to reveal the propagation of the fluctuation. Monitors A/B/C1/D1 are set along the circumference of impeller with a radius of 70.5 mm (the outer diameter of the impeller is 70 mm) to analyze the circumferential variation of pressure. Monitors C1/C2/C3 and D1/D2/D3 are respectively set in the 270◦ and 360◦ directions of the coordinate system. The distance between C1, C2 and C3 is 14.8 mm and the distance between D1, D2 and D3 is 20.8 mm. These monitors are chosen to analyze the time history of static pressure in radial direction between impeller and volute. In order to analyze the influence of volute tongue radius on the surrounding flow, monitors F and E are set at the exit of the impeller close to the volute tongue and the pressure fluctuation data at monitor A is also used for analysis. On the surface of the volute tongue, the apex H and two points G and I where the curved section connects with the straight sections tangentially are selected to analyze the pressure fluctuation of the volute tongue region.

**Figure 18.** Positions for the static pressure monitors within the volute.

#### 3.5.1. Static Pressure Fluctuation along the Circumference of the Impeller Outlet

Figure 19 shows the time history of pressure at monitors A/B/C1/D1 along the circumference of the outlet of the impeller. It is seen that for the *r* = 4 mm model, the pressure exhibits a remarkable fluctuating mode at all four monitors; three periods with large amplitudes are observed corresponding to the three revolutions of the impeller. The remarkable fluctuation shows different patterns at the various monitors. However, for the *r* = 9 mm, 11 mm and 13 mm models, the fluctuating amplitude is greatly reduced for all monitors; the pressure actually fluctuates in a quasi-periodic mode with the number of period during one revolution of impeller equals to the number of blades. The fluctuation at the four monitors shows obvious differences; the amplitude is the maximum for monitor-A and minimum for monitor-B and is mild for monitor-C and monitor-D. This indicates that the unsteadiness of pressure is greatly weakened as the monitor is away from the volute tongue, while the weakening is not proportional to the distance from the monitor to volute tongue but also depends on the interaction of impeller and volute.

**Figure 19.** Static pressure fluctuation at monitors A/B/C1/D1.

The time-averaged value of pressure is also affected the volute tongue radius. For monitor-A, the time-averaged pressure is relatively higher for the *r* = 4 mm and 9 mm models and lower for other two but the difference is minor. The time-averaged value for the *r* = 4 mm model is always the lowest for other three monitors. Although the volute is modified only around the tongue, it does have influence on the flow in the volute far away from the tongue through an indirect way; as the volute tongue radius decreases, the time-averaged pressure at the outlet of the impeller close to the volute tongue is increased, while decreased to a certain extent at other three monitors away from the volute tongue.

Considering the volute geometry exemplified in Figure 18 and the URANS data of radial velocity in Figure 15, it can be found that the radial velocity is near zero in the blade passages corresponding to monitor-A (around 52◦ ), while it increases to around 3 m/s at monitor-B (around 140◦ ) which results in the decreasing pressure from monitor-A to monitor-B. However, the magnitude of radial velocity in the blade passages further increases to 4.5 m/s at monitor-C1 (around 230◦ ) and the pressure is also increased; this may be caused by the significant increase in the distance between impeller and volute. From point C1 to point D1 (around 325◦ ), although the distance is increased to a certain extent than

that at point C1, the radial velocity in the blade passages changes from 4.5 m/s to 8 m/s, so the pressure at point D1 may be more affected by the radial velocity which reduces the local pressure.

For the pressure at the present monitors, the pressure obtained by RANS is normally far deviated from the time-averaged URANS result. The pressure computed by RANS is normally over-predicted for monitor-A, -B and -C1 but is under-predicted by monitor-D and the relative difference can be up to about 30%. For the volute with small tongue radius such as *r* = 4 mm, the minimum relative difference is found at the monitor near the volute tongue, while for other models well consistency is found at monitor-C1 which is the furthest away from the volute tongue.

#### 3.5.2. Static Pressure Fluctuation along the Radial Direction in the Volute

The comparison of pressure fluctuation monitored at C1/C2/C3 reflects the influence of volute tongue radius on flow far away from it, as given in Figure 20. It is seen that generally, the pressure presents strong fluctuation at all three monitors for the *r* = 4 mm model and quasi-periodic fluctuation is clearly seen corresponding to the three revolutions of the impeller. The fluctuation reduces significantly for monitors of other three models. This indicates that although only the volute tongue is modified, it induces the variation of local flow with would further affect the flow in the whole volute in an indirect way.

**Figure 20.** Static pressure fluctuation at monitors C1/C2/C3.

The time-averaged pressure does not vary a lot for different models but is greatly dependent on the radial position of the monitor. For each model, the time-averaged pressure is the lowest for monitor-C1 which locates near the volute outlet and increases as the flow moves outwards toward the volute where the pressure is recovered from the kinetic energy of fluid. For the *r* = 9 mm, 11 mm and 13 mm models, the amplitudes of pressure fluctuation of monitor-C2 and -C3 are similar and are much smaller than that of monitor-C1, which reveals that the fluctuation of near-wall of the volute surface is greatly damped via viscous mechanism by the wall.

By comparing the results obtained by RANS and URANS, it is found that the RANS simulation normally over-predicts the pressure at all monitors. The degree of over-prediction is larger for the model with smaller tongue radius as the flow has stronger unsteadiness and is relatively larger for monitor-C1 since the jet-wake structure for the outflow of impeller is inherently unsteady and could not be revealed by RANS.

The time-history of pressure for flow towards the outlet of the volute is represented by monitors D1/D2/D3 as shown in Figure 21. It is seen that substantial fluctuation of pressure appear for the *r* = 4 mm model, while the fluctuation is greatly attenuated for other models. For each model, the pressure is the highest for the monitor close to the volute surface and lowest close to the impeller, the same variation tendency as the monitors C1/C2/C3. The amplitude of fluctuation is also the largest for monitor-D1 where the unsteady flow out of the impeller dominates. As the volute tongue radius increases, the fluctuating amplitude does not vary much at all in monitors for the *r* = 9 mm, 11 mm and 13 mm models. For the data obtained by monitors D1/D2/D3, it is different from those of C1/C2/C3 that the steady-state RANS simulations always under-predict the time-averaged pressure especially for monitor-D1 which is the closest to the outlet of impeller.

**Figure 21.** Static pressure fluctuation at monitors D1/D2/D3.

3.5.3. Static Pressure Fluctuation Near the Volute Tongue

Since the flow moving out of the impeller will directly impinge on the volute tongue in the region around monitor-A, here we present and discuss the pressure fluctuation at monitors A/F/E in Figure 22 which are located at the outlet of the impeller and around the volute tongue. The time-averaged pressure is the lowest at monitor-A and highest at monitor-E for the *r* = 4 mm model, while for other models the value can be higher for monitor-E or -F depending on the volute radius and the difference is relatively smaller. This is partially attributed to the modification of volute tongue geometry that as its radius decreases, the flow patterns at monitor-E approach to that of monitor-F since monitor-E is also blocked by the volute tongue. Consequently, the position of the highest pressure at the outlet of the impeller near the volute tongue also changes accordingly, that is, distance between the circumferential position at the outlet of impeller with highest pressure and the volute tongue is almost fixed. Although there is a certain degree of difference between RANS and URANS results, we find RANS usually over-predict the pressure at monitor-A and -F and the magnitude of over-prediction gets remarkable as the volute tongue radius increases. For the flow at monitor-E, the predicted pressure of RANS is higher than that of the time-averaged URANS data for the *r* = 4 mm model and lower for *r* = 11 mm and 13 mm models.

**Figure 22.** Static pressure fluctuation at monitors A/E/F.

The monitors G/H/I are positioned right at the surface of volute tongue in which monitor-G and monitor-H are closer to the impeller that the flow may impinge on them, while monitor-I is at the outlet section of the volute and is free of flow impingement. Figure 23 shows the time history of pressure at the three monitors. The pressure fluctuation is still notable for the *r* = 4 mm model, while the amplitude is negligible for the other three models compared with data obtained at other monitors. The time-average pressure for monitor-H at the curved section of the volute tongue has the largest value and is lowest for monitor-G inside of the volute tongue. This indicates that the outflow of the impeller generally impinges on the position of monitor-H but moves partially in a tangential way over the volute surface at monitor-G in the form of recirculating flow. The time-averaged pressure at a certain monitor varies with the volute tongue radius. For all monitors G/H/I, the pressure is the lowest for the *r* = 4 mm model and highest for the *r* = 13 mm model, indicating that the increasing of volute

tongue radius will increase the local pressure at the center of the tongue region; however, the pressure in the volute tongue region close to the volute exit will be reduced to a certain extent.

**Figure 23.** Static pressure fluctuation at monitors G/H/I.

The steady-state RANS simulation gives the worst prediction of pressure especially for monitor-H in the curved section probably due to the inaccurate solution of local flow velocity. It shows that for the pressure distribution at the volute tongue, the high-pressure concentration area predicted by RANS is more biased toward the exit of the volute. For the same model, the RANS results at the monitor-I do not show a pressure variation trend as the URANS results, while it is the same for monitor-G and monitor-H.

#### *3.6. Aeroacoustic Characteristics*

In this work, the internal flow of the Sirocco fans is numerical investigated using the RANS and URANS approaches considering the efficiency and feasibility of computational resources. The results reflect the general characteristics of the temporal variation and spatial distributions of certain quantities but is insufficient in capturing the subtle pressure fluctuation, thus the aeroacoustic characteristics cannot be numerically computed. In this section, we present the experimental results regarding the aeroacoustic characteristics of the various models.

The influence of volute tongue radius of the far-field noise characteristics are analyzed in this section. The schematic of the arrangements of microphones is shown in Figure 24; four microphones A/B/C/D are positioned at the inlet of the fan and one microphone-E is positioned at the outlet of the fan. The sound pressure level (SPL) measured at each microphone is listed in Table 5 for the various volute models. For the microphones A/B/C/D at the outlet of the fan, it is seen in the table that the SPL measured at microphones A/B/C/D increases significantly as the volute tongue radius decreases especially for the *r* = 4 mm model. The SPL of the *r* = 9 mm, 11 mm and 13 mm models is of the same level with maximum relative difference less than 2%. The SPL measured at microphone-E at the inlet of the fan exhibits smaller difference for the several models. By comparing the SPL at the several microphones, it can be found that the SPL measured by microphone B/C is relatively higher than that of A/D because the flow at the outlet of the volute has a larger velocity in the direction of B/C which will produce a stronger SPL.

**Figure 24.** Positions of microphones in the aeroacoustic test. (**a**) Top view of the positions of the probes; (**b**) side view of the positions of the probes.


**Table 5.** The sound pressure level monitored by different microphones (dB).

Figure 25 shows the spectrum of 1/3 octave band of noise obtained at monitor-B at the outlet of volute and monitor-E at the inlet of volute. Considering that the blade passage frequency (BPF) of the impeller is 820 Hz, it is seen in the figure that the maximum SPL of the two monitors is around the frequency of BPF for all models, thus it contributes the most to the overall SPL of the fan. The volute model with smaller radius has a higher magnitude of SPL roughly for the low frequency regime, while in the high frequency regime the volute with larger radius has a higher SPL. The SPL increases significantly in most of the low frequency regime for volutes with small tongue radius, while the reduction in the high frequency regime is relatively small. As the volute tongue radius varies within a certain range, such as *r* = 9 mm to *r* = 13 mm, the SPL in most frequency regimes is relatively less affected by the tongue radius.

**Figure 25.** Noise spectrum: (**a**) Spectrum of 1/3 octave band of microphone-B; (**b**) Spectrum of 1/3 octave band of microphone-E.

### **4. Conclusions**

This work presents an experimental and numerical investigation on the influence of volute tongue radius on the aerodynamic and aeroacoustic characteristics of a Sirocco fan. The simulations are performed by both RANS and URANS approaches to comparably explore their capability, accuracy and reliability for such configurations. The characteristics of internal flow under the effect of different volute tongue radius are presented in terms of the aerodynamic quantities, the flow within the impeller and near-wall flow on the volute surface, the pressure fluctuation within the volute and the aeroacoustic characteristics. The experimental and numerical results have the following conclusions:


It should be emphasized here that the findings concluded in this work is considered only suitable for fan models of similar geometry, especially the small-size Sirocco fans and centrifugal fans under certain flow conditions. The effect of volute tongue radius on other types of fans or fluid machines needs to be investigated in the following works. The comparison of RANS and URANS approaches is also limited to the *k*-ε turbulence model and the conclusions may be different both in variation trend and magnitude as other turbulence models are employed.

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

**Funding:** This work was supported by National Natural Science Foundation of China (51706205), Zhejiang Public Welfare Project (LGG20E060001) and Zhejiang Province Science and Technology Plan Project (2020C04011).

**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* **Large Eddy Simulation of Film Cooling with Triple Holes: Injectant Behavior and Adiabatic Film-Cooling E**ff**ectiveness**

#### **Seung Il Baek and Joon Ahn \***

School of Mechanical Engineering, Kookmin University, 77 Jeongneung-ro, Seongbuk-gu, Seoul 02707, Korea; greenjet50@gmail.com

**\*** Correspondence: jahn@kookmin.ac.kr

Received: 19 October 2020; Accepted: 9 November 2020; Published: 11 November 2020

**Abstract:** This study investigated the effect of adding two sister holes placed downstream the main hole on film cooling by employing large eddy simulation. Here, film-cooling flow fields from a triple-hole system inclined by 35◦ to a flat plate at blowing ratios of M = 0.5 and unity were simulated. Each sister hole supplies a cooling fluid at a flow rate that is a quarter of that for the main hole. The simulations were conducted using the Smagorinsky–Lilly model as the subgrid-scale model, and the results were compared with those for a single-hole system for the same amount of total cooling air and same cross-sectional area of the holes. Relative to the single-hole system, the spanwise-averaged film-cooling effectiveness in the triple-hole configuration at M = 1.0 increased by as much as 345%. The subsequent proper orthogonal decomposition analysis showed that the kinetic energy of a counter-rotating vortex pair in the triple-hole system dropped by 30–40% relative to that of the single-hole system. This indicates that the additional sister holes promoted the suppression of the mixing of the coolant jet with the mainstream flow, thus keeping the temperature of the latter low. Cross-sectional views of the root-mean-square temperature contours were also analyzed; with the results confirming that the effect of the sister holes on the jet trajectory greatly contributes in promoting film-cooling effectiveness as compared to the effect of the reduced mixing.

**Keywords:** film cooling; large eddy simulation; triple holes; blowing ratio; adiabatic film-cooling effectiveness; proper orthogonal decomposition

#### **1. Introduction**

The principle of idealized Brayton cycle states that gas turbine efficiency can be improved by increasing the inlet temperature of the turbine [1]. Here, to avoid the excessive thermal stress build-up on the turbine blades, the surface temperature of the blades must remain below a tolerable limit. Various cooling methods can be employed, among which the film-cooling technique has been widely used. In this technique, the cooling air which is bled from the compressor is injected through holes on the turbine blade surface—this reduces the surface temperature and protects the surface from the main flow. For the holes, a cylindrical shape is often appropriate because of its simplicity and ease of manufacture.

However, substantial investigations of the physics of film cooling in such a cylindrical hole system have revealed that a single cylindrical hole can be relatively vulnerable to the generation of the jet liftoff [2] and that a specific hole arrangement could produce film cooling that provides better thermal protection than single cylindrical holes do. Practically, specific hole arrangements make the intended vortex interactions between film cooling jets to increase film-cooling effectiveness.

Film cooling in turbine engines can be effectively studied at a reasonable cost using computational fluid dynamics (CFD). High turbulence is generated in the film-cooling flow field, after which a model is constructed for the turbulence. Another option is the well-known large eddy simulation (LES), which can predict the mixing between the main flow and coolant jet better than Reynolds-averaged Navier–Stokes Simulation (RANS) [3–5], although such requires much higher computational costs. Comparatively speaking, RANS cannot accurately predict the complex flow structures induced in film-cooling flow fields as all of the turbulent fluctuations are ensemble-averaged, whereas LES resolves large-scale eddies directly in the turbulent flow, resulting in more accurate predictions of the complex flow [6,7].

Heidmann et al. [8] showed in a numerical investigation that relative to a one-hole system, three holes (consisting of the main hole and two side holes) can increase film-cooling effectiveness, as well as decrease the heat transfer coefficients at a blowing ratio of unity, using RANS. They also reported that the fabrication of cylindrical triple holes is economically more viable than that of other shapes. Meanwhile, Javadi et al. [9] constructed square sister holes to control vortex interactions between the jet and crossflow at M = 0.5 by employing RANS. They found that two sister holes with a square cross-section reduce the mixing between the coolant and crossflow, thus improving film-cooling effectiveness. Ely et al. [10] also used RANS to investigate an increase in film-cooling effectiveness with two cylindrical sister holes. Using blowing ratios of 0.2, 0.5, 1.0, and 1.5, they demonstrated that definite improvements in film-cooling effectiveness occur at high blowing ratios.

Furthermore, Wu et al. [11] conducted both experimental and numerical studies of film cooling given three sister-hole configurations and confirmed an improvement in the performance as above. Here, they compared the experimental and numerical data only in terms of spatially averaged film-cooling effectiveness, which showed a difference of more than 30% under blowing ratios of unity and 1.5.

Farhadi-Azar et al. [12] introduced LES to study the effects of the triple jets on film cooling. On the basis of their demonstration, triple jets from the main hole along with two small sister holes increased film-cooling effectiveness because of weaker counter-rotating vortex pairs in the jet from the main hole, which led to less mixing between the coolant jets and main flow. However, several limitations could be found in their study. For instance, they employed an injection angle (α) of 90◦ on their hole, which is very uncommon for a gas turbine design because of the strong jet liftoff. Aside from this, the plenum was not included in the simulation, and a flow with a uniform velocity was applied to the hole inlet, which is unrealistic in gas turbine scenarios.

Considering the above, this paper aims to investigate the effect of a triple-hole system on film-cooling performance on a flat plate at blowing ratios of 0.5 and 1.0. The simulation was conducted using LES under a more realistic flow condition and at an injection angle (α) of 35◦ for the hole. The plenum is included, and the velocity in the hole is nonuniform with both the high-momentum region generated by the jetting effect and low-momentum region.

The LES results are used to understand the effect of the flow generated in the triple hole on the trajectory of the film-cooling jet and mixing with the mainstream flow. Vortex structures that affect the mixing are identified through time-averaged velocity fields and instantaneous flow fields, whereas the contribution of these vortices is determined through a proper orthogonal decomposition (POD) analysis. The film-cooling performance, degree of mixing with the main flow, and trajectory of the film-cooling jet are also observed through the time-averaged temperature fields. On the basis of the analysis results, the effect and mechanism of the triple hole on the film-cooling performance are identified.

#### **2. Geometry and Boundary Conditions**

The geometry of the single-hole systemin the test casewas taken from Seo et al. [13], whose experimental apparatus consisted of a row of five cooling holes. Herein, a single-hole configuration was adopted to simplify the computational domain, and a periodic boundary condition was applied on the mesh sides to reduce the computational cost. The computational domains of the single- and triple-hole systems are delineated by the orange dashed lines in Figure 1, whereas a schematic of the cooling holes in each case is shown in Figure 2. In the single-hole geometry (Figure 2a), the hole diameter *D* was 25 mm. The triple holes comprised two sister holes which are placed slightly downstream of the primary hole (Figure 2b) as described in Ely et al. [10]. In the triple-hole geometry, the hole diameters were calculated to match their total cross-sectional areas to the cross-sectional area of the single-hole. Thus, the primary hole diameter *D*′ was 20.4124 mm, and the diameter of each sister hole was *D*′ /2 = 10.2062 mm. The hole length-to-diameter ratio (*L*/*D*), injection angle (α), and pitch-to-hole diameter *(P*/*D*) were set to 1.6, 35◦ , and 3.0, respectively. In addition, the holes were cylindrical, with no compound angles. ′ ′ α ′ α

′

**Figure 1.** Computational domains of the single- and triple-hole configurations (orange dashed lines).

**Figure 2.** Schematics of the single hole and triple holes cooling hole configurations.

The boundary conditions in each computational domain are specified in Table 1. The turbulence intensity of the mainstream flow was 0.2% at the main inlet (as reported experimentally [13]), and the main flow velocity was 10 m/s. The temperatures of the main flow and coolant at the inlets were 313 and 293 K, respectively. Moreover, the boundary layer thickness at the cooling hole center is about *D*.

**Table 1.** Boundary conditions of the computational fluid dynamics (CFD) domain.


As is well-known, injecting the film-cooling jet induces a counter-rotating vortex pair (CRVP) in the system (see Figure 3a). This vortex pair is sourced from the vortices generated in the boundary layer of the coolant tube wall and around the leading edge of the hole [12,14,15]. Interactions between the jet and main flow cause the injected coolant to split into two directions, whereas the vortices generated in the coolant are stretched and turned. As the CRVP vortices draw closer together, their mutual induction increases, promoting coolant jet liftoff. The consequent strong entrainment of the hot main flow under the jet increases the adiabatic temperature on the test plate and decreases film-cooling effectiveness. Moreover, when two sister holes are placed slightly downstream of the primary hole (Figure 2b), their jets generate an anti-CRVP, as illustrated in Figure 3b. Therefore, film-cooling effectiveness is increased by establishing vortex interactions between the jets and by reducing the intensity of the main CRVP [9–12].

**Figure 3.** Cross-sectional schematic of a counter rotating vortex pair generated in the single- and triple-hole configurations.

#### **3. Validation of Numerical Methods**

The CFD calculations were performed in ANSYS Fluent v.19.1 [16], and the meshes were generated in Pointwise v.18.1 [17]. The numerical simulations were executed in the LES using the Smagorinsky–Lilly model as the subgrid-scale model. The time step was set to 6.25 <sup>×</sup> <sup>10</sup>−<sup>6</sup> s. Assuming the mainstream convected the hole diameter after 400 time steps, the simulations were carried out with the computational time step <sup>∆</sup>*<sup>t</sup>* = 0.0025*D*/*U*<sup>∞</sup> [18,19]. Once the system had reached a statistical steady state, the statistics of the solution was collected in multiples of the period. In each time step, approximately 10 subiterations were executed to obtain well-resolved data [20]. The LES runtime using 18 cores of an Intel Xeon Gold 6148 processor was 1–2 months in each case.

The fluid was assumed as Newtonian and incompressible and had temperature-dependent variable properties. The compressibility effect was negligible because the main flow velocity was 10 m/s (Mach 0.029) and the jet injection velocities were 5 and 10 m/s at *M* = 0.5 and *M* = 1.0, respectively; that is, both velocities were far below Mach 0.3 [21]. The governing equations of the CFD simulations are the continuity and momentum equations. The filtered Navier–Stokes equations in the LES approach are given by [22]

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_i} (\rho \overline{u}\_i) = 0 \tag{1}$$

and

*τ*

$$\frac{\partial(\rho \overline{u})}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (\rho \overline{u}\_i \overline{u}\_j) = \frac{\partial}{\partial \mathbf{x}\_j} \Big[ \mu \Big( \frac{\partial \overline{u}\_i}{\partial \mathbf{x}\_j} + \frac{\partial \overline{u}\_j}{\partial \mathbf{x}\_i} \Big) - \frac{2}{3} \mu \frac{\partial \overline{u}\_l}{\partial \mathbf{x}\_l} \delta\_{ij} \Big] - \frac{d \overline{p}}{d \mathbf{x}} + \frac{\partial \tau\_{ij}}{\partial \mathbf{x}\_j} \tag{2}$$

where τ*ij*, the subgrid-scale turbulent stress, needs modeling using the Boussinesq hypothesis like RANS models as

∆ = 0.0025/<sup>ஶ</sup>

$$
\pi\_{ij} - \frac{1}{3}\pi\_{kk}\delta\_{ij} = -\mu\_l \left(\frac{\partial \overline{u\_l}}{\partial \mathbf{x}\_j} + \frac{\partial \overline{u\_j}}{\partial \mathbf{x}\_l}\right) \tag{3}
$$

where µ*<sup>t</sup>* is the turbulent viscosity on the subgrid scale. In the Smagorinsky–Lilly model, the turbulent viscosity is modeled as [18]

$$\mu\_l = \rho L\_s^{-2} \sqrt{\frac{1}{2} \left( \frac{\partial \overline{u\_l}}{\partial \mathbf{x}\_j} + \frac{\partial \overline{u\_j}}{\partial \mathbf{x}\_i} \right) \left( \frac{\partial \overline{u\_l}}{\partial \mathbf{x}\_j} + \frac{\partial \overline{u\_j}}{\partial \mathbf{x}\_i} \right)} \tag{4}$$

(ത) <sup>+</sup> ൫തത൯ = ቈ ቆഥ<sup>ప</sup> + ഥ<sup>ఫ</sup> ቇ − <sup>2</sup> 3 ഥ − ̅ <sup>+</sup> Figures 4 and 5 show the meshes of the single- and triple-hole systems in the *xy* and *yz* planes, respectively. The computational domains of the systems were composed of 2.04 (single-hole) and 2.69

ഥ<sup>ప</sup> + ഥ<sup>ఫ</sup> ቇ

= −<sup>௧</sup> ቆ

 − 1 3 *μ*

*μ*

(triple-hole) million hexahedron cells. Figure 6 shows a close-up of the mesh near the hole, where P2 and P1 are the static pressures at the inlet and outlet of the hole, respectively, as mentioned earlier.

 ቇ ቆഥ<sup>ప</sup> + ഥ<sup>ఫ</sup> ቇ

 ቇ ቆഥ<sup>ప</sup> + ഥ<sup>ఫ</sup> ቇ

<sup>௧</sup> = <sup>௦</sup>

<sup>௧</sup> = <sup>௦</sup>

ଶඨ 1 2 ቆ ഥ<sup>ప</sup> + ഥ<sup>ఫ</sup> 

ଶඨ 1 2 ቆ ഥ<sup>ప</sup> + ഥ<sup>ఫ</sup> 

(**b**) Triple holes

**Figure 4.** CFD meshes of the single- and triple-hole systems in the *xy* plane.

(**a**) Single hole (**b**) Triple holes

**Figure 6.** Close-up of the mesh near the hole.

Figures 7 and 8 show the results of the grid sensitivity tests in each system at a blowing ratio of 0.5, with the grid descriptions given in Tables 2 and 3. The results were obtained using the Smagorinsky–Lilly model in LES, and the tests were performed on five different grids. In the single-hole system, the

effectiveness value along the centerline (*z* = 0, *y* = 0) on the third grid almost matched those of the finer fourth and fifth grids.

**Figure 7.** Grid sensitivity test of the large eddy simulation (LES) calculation in the single-hole system.

**Figure 8.** Grid sensitivity test of the LES calculation in the triple-hole system.

**Table 2.** Specifications of mesh arrangements in the grid sensitivity test of the single-hole system.


**Table 3.** Specifications of mesh arrangements in the grid sensitivity test of the triple-hole system.


Therefore, the adiabatic film-cooling effectiveness in the single-hole system was evaluated on the 2.04 million-cell grid. In the triple-hole system, the third grid with 2.69 million cells again showed almost the same centerline effectiveness as the fourth and fifth grids. Therefore, the 2.69 million-cell grid was selected for modeling the film-cooling effectiveness in the triple-hole system.

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

#### *4.1. Centerline and Spanwise-Averaged Film Cooling E*ff*ectiveness*

Figure 9 shows the centerline and spanwise-averaged effectiveness values in the single- and triple-hole systems at *M* = 0.5 and 1.0. The centerline effectiveness decreased as the main flow traveled downstream. This degradation is attributable to the turbulence generation and rapid mixing between the main flow and jet, which increases the test plate temperature. For the same total cross-sectional hole area and same amount of total cooling air, the triple holes improved the centerline and spanwise-averaged effectiveness values, respectively, by approximately 36% and 45% over the single-hole configuration at *M* = 0.5 and approximately 250% and 345% at *M* = 1.0. In Figure 9a–d, the LES-predicted centerline film-cooling effectiveness in the single-hole system showed a good match to the experimental data by Seo et al. [13], whereas the LES model under-predicted the triple-hole, spanwise-averaged effectiveness experimental data by Cao et al. [22] because the data were measured for *L*/*D* = 10. Thus, compared to lower *L*/*D* cases, higher *L*/*D* values of the hole increase the adiabatic film-cooling effectiveness under the same operating condition as the coolant velocity at the hole exit becomes more uniform [13].

**Figure 9.** Film cooling effectiveness values of the single- and triple-hole systems.

The ratios of the spanwise-averaged effectiveness values in the triple-hole system to those in the single-hole system as functions of *x*/*D* for *M* = 0.5 and 1.0 are plotted in Figure 10. Note that increasing the blowing ratio improves the ratio of the spanwise-averaged effectiveness. This indicates that under high blowing ratios, configuring the main hole and two sister holes is more effective for film cooling than using an ordinary single hole mainly because the jet trajectory in the single-hole system is mostly governed by the jet liftoff at a high blowing ratio of *M* = 1.0. The result is very low adiabatic film-cooling effectiveness, as indicated in Figure 9c,d, and very high ratios of the spanwise-averaged effectiveness values, as shown in Figure 10. The maximum value of the ratio at *M* = 1.0 is observed around at *x*/*D* = 4 since the minimum value of the spanwise-averaged effectiveness in the single-hole system is shown around at *x*/*D* = 4 due to the jet lift off, as shown in Figure 9d.

**Figure 10.** Ratios of spanwise-averaged effectiveness in the triple hole system to those in the single-hole system.

#### *4.2. Contours of Film Cooling E*ff*ectiveness and Dimensionless Temperature*

Figure 11 presents top views of the time-averaged film-cooling effectiveness contours on the test plate in each hole system. The additional sister holes largely increased the amount of coolant in contact with the test plate, which boosts film-cooling effectiveness. Moreover, the spread of the coolant in the spanwise direction in the triple-hole system was much larger than that in the single-hole configuration as the main CRVP weakens. Furthermore, lower film-cooling effectiveness was predicted at *M* of unity than at *M* = 0.5 because of strong jet liftoff.

**Figure 11.** Time averaged film cooling effectiveness on the test plate.

The cross-sectional views of the dimensionless mean temperature contours at *x*/*D* = 2, 5, 10, 15, and 20 in the two systems are illustrated in Figure 12. Initially, the dimensionless temperature near the test plate was higher in the triple-hole system than in the single-hole configuration, which increased film-cooling effectiveness on the test plate. As the main flow went downstream, the dimensionless temperatures near the test plate in the triple-hole system became similar to those in the single-hole system, leading to almost identical film-cooling effectiveness in the range of *x*/*D* > 15 (Figures 8 and 9) because of intensive mixing between the mainstream flow and coolant jet.

**Figure 12.** Cross-sectional contours of the dimensionless mean temperatures.

Figure 13 shows the central cross-sectional contours of the dimensionless mean temperatures in the *z* = 0 plane for both hole systems. The entrainment of the main flow under the cooling jet in the triple-hole system was noticeably weaker than that in the single-hole system. Many researchers have argued that the jet lifted off because of the high momentum of the coolant and reattached downstream to the plate at a high blowing ratio [2,22,23]. However, as illustrated in Figure 13b, the dispersion of the coolant jet increases without reattachment with the test plate as the flow goes downstream after the jet liftoff, followed by the entrainment of the main flow under the jet. In addition, the dimensionless temperature near the test plate in the near-hole region in the triple-hole system was much higher than that in the single-hole configuration as vortex interactions between the jets were established, thereby reducing the intensity of the main CRVP at *M* = 0.5 and 1.0 and increasing film-cooling effectiveness on the test plate.

**Figure 13.** Central cross-sectional contours of the dimensionless mean temperatures in the *z* = 0 plane.

Figure 14 illustrates the streamlines of cooling jets to a cross-sectional plane at *x*/*D* = 5 and film-cooling effectiveness on the test plate in the single- and triple-hole systems at *M* = 0.5 and unity. Unlike in a single hole, more streamlines of the coolant head to the region near the test plate in the *x*/*D* = 5 plane appeared in the triple-hole system, with these streamlines more dispersed, which resulted in higher adiabatic film-cooling effectiveness on the plate given that most streamlines from the sister holes have converged to the centerline. At *M* = 1.0, the strong jet liftoff caused the streamlines to cover smaller regions laterally in the *x*/*D* = 5 plane even in the triple-hole system.

**Figure 14.** The streamlines of cooling jets to cross-sectional plane at *x*/*D* = 5.

#### *4.3. Contours of x-Vorticity and Mean Velocity Vectors*

Figure 15 shows the cross-sectional contours of the mean *x*-vorticity at the *x*/*D* = 2 and *x*/*D* = 5 planes in the single- and triple-hole systems at *M* = 0.5 and unity. The mean velocity vectors are superimposed on the contours. As shown in Figure 15b,d of the triple-hole system, distinct anti-CRVPs were found from the sister holes around the *x*/*D* = ±0.5 and *y*/*D* = 0.2 regions in the *x*/*D* = 2 plane. The anti-CRVPs weakened the main CRVP and reduced the main jet liftoff, thereby reducing the entrainment of the hot mainstream gas under the coolant and increasing the dimensionless mean temperature and film-cooling effectiveness, as illustrated in Figures 11 and 13. However, the intensity of the anti-CRVPs became very small in the *x*/*D* = 5 plane because the anti-CRVP from the sister holes was not as strong as the main CRVP. In the triple-hole system, the upward velocity vectors around the line of *z*/*D* = 0 were less concentrated than those in the single-hole system. Thus, coolant liftoff was weakened by the generated anti-CRVPs. At *M* = 1.0, the positions of the core for the CRVPs were higher than those at *M* = 0.5 because of the strong jet liftoff generated.

**Figure 15.** Cross-sectional contours of the *x*-vorticity and velocity vectors.

Figure 16 shows the cross-sectional velocity vectors of POD Mode 1 at *x*/*D* = 2 in the two-hole systems at M = 0.5 and unity. The POD was introduced by Dr. John Lumley in 1967 in order to decompose a random velocity vector field of turbulence into a set of portions of the total kinetic energy [24,25]. On the basis of the figure, CRVP was most distinct in the velocity vector field at *x*/*D* = 2. Figure 17 illustrates the portion of the total kinetic energy at *x*/*D* = 2 for Modes 1, 2, 3, and 4 in the single- and triple-hole systems at *M* = 0.5 and *M* = 1.0, which showed Mode 1 as the most dominant. At *M* = 0.5, POD Mode 1 contributed to about 60% and 36% of the total kinetic energy in the single- and triple-hole systems, respectively, which means that CRVP was dominant for both, although relative to the former, the portion of the kinetic energy in the latter system was largely decreased by about 40% because of a weakened main CRVP by anti-CRVPs from the two sister holes. Similarly, at *M* = 1.0, POD Mode 1 contributed to about 71% and 49% of the total kinetic energy in the single- and triple-hole systems, respectively, indicating that the portion of the kinetic energy was decreased by nearly 30%. Indeed, it could be inferred that compared to that in the single-hole system, the kinetic energy of the CRVP in the triple-hole system was reduced by about 30–40%.

**Figure 17.** Portions of the total kinetic energy for Mode 1, 2, 3, and 4 in the plane of *x*/*D* = 2.

#### *4.4. Root-Mean-Square Temperature Contours*

Figure 18 illustrates the cross-sectional views of the root-mean-square (RMS) temperature contours at *x*/*D* = 2, 5, 10, 15, and 20 for the mixing region between the mainstream flow and coolant in the single- and triple-hole systems. The mixing area was more expanded in the triple-hole system in the spanwise direction, especially in the near-hole region as *x*/*D* = 2, because the weakened main CRVP allowed wider spreading and closer contact of the coolant jet with the test plate. Although no significant difference in the mixing intensity was observed in both hole systems, the *Trms* in the *x*/*D* = 2 plane reached maximum values of 6.25, 6.24, 6.08, and 5.19, as shown in Figure 18a–d, respectively, which indicates a similar but slightly decreased mixing intensity in the triple-hole system relative to the single-hole system. This reduced mixing intensity is beneficial to film cooling because increased mixing between the coolant and mainstream flow consequently decreases film-cooling effectiveness [22]. Nevertheless, it could be stated that in the triple-hole system, the effect of the sister holes on the jet trajectory contributes greatly to improve film-cooling effectiveness as compared to the effect of less mixing between the main flow and coolant.

**Figure 18.** Cross-sectional contours of the *Trms.*

#### *4.5. Q-Criterion Iso-Surfaces*

In the Q-criterion, Q is defined as the difference between the magnitudes of the rotation rate and strain rate, representing the dominance of the former over the latter [26–28]. The instantaneous contours of a Q-criterion with iso-surfaces of level 80,000 at *M* = 0.5 and unity are displayed in Figure 19. Q-criterion contours in the single-hole system clearly showed various vortex structures, such as the CRVP, hairpin and jet shear layer vortices, and horseshoe vortex, whereas those in the triple-hole configuration had smaller hairpin and jet shear layer vortices and were not very distinguishable. These smaller and indistinct hairpin vortices represent the weakened CRVP as the hairpin vortex is closely related to the generation of the CRVP [22]. A weakened CRVP positively affects the cooling performance on the plate because it increases film-cooling effectiveness on the test plate.

**Figure 19.** Q-criterion contours with iso-surfaces of level 80,000.

#### *4.6. Cross-Sectional Contours of the Mean Velocity Magnitude in the Hole*

− Figure 20 shows the cross-sectional contours of the mean velocity magnitude in the hole measured around the hole exit (*y*/*D* = −0.2). As mentioned earlier, the ratio of the delivery tube length to the hole diameter in the single-hole system was 1.6, whereas it was 1.8 and 3.6 for the main hole and sister holes in the triple-hole system, respectively, because the triple holes had the same total cross-sectional hole area as the single hole. Therefore, the mean velocity magnitude in the main hole is not uniform but has high- and low-momentum regions, whereas that in the sister holes shows more developed coolant flow. The cooling jet having a uniform velocity profile is advantageous to film-cooling performance because of less jet liftoff. Thus, it could be stated that relative to the single-hole system, the sister holes improve film-cooling effectiveness by the generation of the anti-CRVP, as well as higher *L*/*D* ratios.

−

−

− **Figure 20.** Cross-sectional contours of the mean velocity magnitude in the hole at *y*/*D* = −0.2.

#### **5. Conclusions**

Film-cooling simulations of the single- and triple-hole systems on a flat plate were carried out in this paper at blowing ratios of *M* = 0.5 and unity by employing LES. On the basis of these simulations, when two sister holes which were placed slightly downstream the main hole, the coolant jets from the former generated an anti-CRVP, which increased film-cooling effectiveness. Thus, for the same total cross-sectional hole area and same amount of total cooling air, the triple-hole system increased the centerline and spanwise-averaged effectiveness values by approximately 36% and 45% at M = 0.5 and by approximately 250% and 345% at *M* = 1.0 as compared to that done by the single-hole configuration. Moreover, the ratios of the spanwise-averaged effectiveness values in the triple-hole system indicate that it is the more effective system for film cooling under a high blowing ratio.

Meanwhile, the time-averaged film-cooling effectiveness contours on the test plate and the cross-sectional views of the dimensionless mean temperature contours in the two hole systems showed that additional sister holes can increase the amount of coolant contacting with the test plate and boost film-cooling effectiveness. Specifically, the mean dimensionless temperature contours in the *z*/*D* = 0 plane showed that in the single-hole system at *M* of unity, as the cooling jet goes downstream, the jet gets dispersed without reattachment with the test plate, and the entrainment of the main flow under the jet is generated.

Furthermore, the cross-sectional mean x-vorticity contours in the triple-hole system showed that the generated anti-CRVPs weaken the main CRVP. The results of the POD analysis confirmed that the kinetic energy of the CRVP in the triple-hole system was reduced by about 30–40% relative to that in the single-hole system. Streamlines of the cooling jets to the cross-sectional plane at *x*/*D* = 5 also showed that most streamlines from the sister holes converge to the centerline. As with the Q-criterion contours in the triple-hole system, the hairpin and jet shear layer vortices were smaller and indistinguishable, unlike in the single-hole system.

Lastly, on the basis of the cross-sectional views of the RMS temperature contours, the effect of the sister holes on the cooling jet trajectory contributes significantly to the increase in film-cooling effectiveness, as compared to the effect of less mixing. Accordingly, the sister holes improve film-cooling effectiveness using higher *L*/*D* ratios.

**Author Contributions:** S.I.B. carried out the simulations, analyzed the CFD data and wrote the paper. J.A. supervised the research, analyzed the CFD data and wrote the paper. 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 declared no conflict of interest.

### **Nomenclature**


Θ dimensionless temperature =

Subscripts


#### **References**

1. Moran, M.; Shapiro, H.; Boettner, D.; Bailey, M. *Fundamentals of Engineering Thermodynamics*, 8th ed.; Wiley: Hoboken, NJ, USA, 2014; p. 532.

(*TG*−*T*) *TG*− *T<sup>C</sup>*


**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* **Influence of Tip Clearance on Flow Characteristics of Axial Compressor**

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

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

**\*** Correspondence: byang0626@sjtu.edu.cn; Tel.: +86-21-3420-6871

Received: 9 September 2020; Accepted: 10 November 2020; Published: 11 November 2020 -

**Abstract:** This paper studies the influence of tip clearance on the flow characteristics related to the performance. Based on full-passage numerical simulation with experimental validation, several clearance models are established and the performance curves are obtained. It is found that there exists an optimum clearance for the stable working range. By analyzing the flow field in tip region, the role of the tip leakage flow is illustrated. In the zero-clearance model, the separation and blockage along the suction side is the main reason for rotating stall. As the tip clearance is increased to the optimum value, the separation is suppressed by the tip leakage flow. However, with the continuing increasing of the tip clearance, the scale and strength of the tip clearance vortex is increased correspondingly. When the tip clearance is larger than the optimum value, the tip clearance vortex gradually dominates the flow field in the tip region, which can increase the unsteadiness in the tip region and trigger forward spillage in stall onset.

**Keywords:** axial compressor; tip clearance; flow field; numerical simulation

#### **1. Introduction**

It is the tip region of a rotor that is one of the most unstable regions in an axial compressor. For example, for the tip-critical rotors, the rotating stall was first observed in the tip region. In the tip region, there exists the maximum rotating linear velocity, tip leakage flow, and boundary layer flow, etc. Among these flows, the tip leakage flow, induced by the pressure difference between blade sides, plays the most important role on the flow behavior. The tip leakage flow is closely related to the tip clearance, which is one of the key parameters for the design and manufacture of an axial compressor. As the tip clearance is changed, not only is the strength of the leakage flow changed, but also other flow structure and distribution are influenced greatly, especially under a near-stall condition. Therefore, the performance and stable working range is varied with the tip clearance.

In the last century, many researches focused on the influence of tip clearance in an overall manner. Combined experiment data with potential-vortex model, Lakshminarayana [1] proposed a formula, firstly, to evaluate the stage efficiency in different tip clearance. Based on this formula, the pressure ratio and efficiency can be improved by the reduction of the tip clearance. Inoue et al. [2,3] obtained and described the roll-up tip clearance vortex(TCV) by hotwire. From the vortex model, the flows in the tip region could be divided into the main flow region and tip leakage flow region. Day [4] also proposed the influence of the tip clearance on off-design conditions in axial compressors.

As the experiment method and computational fluid dynamic method have been both highly improved, the flow structures and special phenomenon in the tip region were studied in detail, which included double leakage, vortex breakdown, self-induced flow, and active control technology, etc. For double leakage, it was first noticed by Khalsa et al. [5] in the experiment and simulation results. They pointed out that the relative stagnation pressure is decreased by double leakage; therefore, the loss

and blockage were increased. Double leakage flow was much influenced by clearance size, working condition, and unsteady interaction from the stator, which was illustrated by Mailach et al. [6,7], Sirakov et al. [8], and Du et al. [9]. By the means of Large-eddy simulation (LES) approach, Hah [10] studied the double leakage flow under the large tip clearance model. As double leakage flow occurred in all working conditions with large tip clearance, vortex rope is generated near leading edge, which results in the high mixing loss and non-synchronous vibration.

For the breakdown of tip clearance vortex, Lotzerich [11] and Furukawa et al. [12] both noticed this phenomenon at almost the same time. Later, Furukawa et al. [13,14] summarized the characteristics by the means of experiment and computational fluid dynamics (CFD). When the vortex breakdown occurs, the disappearance of vortex core and the roll-up structure causes the large low-velocity cells. The vortex breakdown always occurred in near-stall condition, because the adverse pressure gradient and the unsteadiness reached the maximum values in the tip region. Yamada et al. [15] presented the unsteady flowing field near the tip region by the CFD method with a single passage grid. It is shown that the TCV is broken down by the effect of the shock wave, and the type of vortex breakdown in near-peak point was different with its in near-stall condition.

From the research above, it is understood that, under near-stall condition, the tip clearance vortex will be highly unsteady and unstable. Therefore, the stall inception was always detected in the tip region. Using a hotwire arranged around the rotor tip of the axial compressor, circumferentially, McDougall et al. [16,17] captured the large propagating wave with small amplitude in near-stall condition, which is called a modal type rotating stall. Later, Day et al. [4,18] captured the short propagating wave with large amplitude, which was a spike type rotating stall. Hoying et al. [19] studied the relation between the rotating stall and tip clearance vortex in the E3 compressor. By using Reynolds averaged Naiver-Stokes simulation (RANS) method, it was found that the tip clearance vortex will move forward to rotor inlet in near-stall condition. Hah et al. [20] and Vo et al. [21] proposed the criterion for the spike initiated rotating stall, which had leading edge spillage and trailing edge backflow. This criterion was already accepted by many researchers. However, many people misunderstood this criterion and considered the tip clearance vortex as the origin of the rotating stall.

By means of unsteady multi-passage simulation, Pullan et al. [22] illustrated two different stall mechanisms. It was pointed out that the stall originated from high incidence of the rotor inlet. Although tip clearance was not necessary for the rotating stall, it does influence the stall characteristics a lot. Later, Smith et al. [23] illustrated the influence of tip clearance on the rotating stall. There existed an optimum tip clearance for stall margin, which can balance the effect from the different stall mechanisms.

However, in their works, the stall was triggered by several re-staggered blades. It was more appropriate to use the non-re-staggered blade model. Although they studied the relationship between stall and tip clearance, they focused on the stall mechanism in the zero clearance model, the evolution of stall mechanism and the effect of the tip clearance vortex in different clearances were not analyzed enough. In addition, the benefit of leakage flow was not mentioned.

This paper studies the influence of tip clearance on flow characteristics related to performance. Based on the full-passages model, steady and unsteady results are obtained by the RANS method. The relationship between the stability of the compressor and the tip clearance is discussed. Then, the flowfield characteristics under different tip clearances are compared, and the effects of the separation and the tip leakage flow on the efficiency and stability are discussed.

#### **2. Experiment and Numerical Setup**

#### *2.1. Description of the Test Rig and the Measurement*

The test was conducted in the 1–1/2 stage axial subsonic compressor in Shanghai Jiao Tong University. The blade number of inlet guide vane (IGV) is 32, and 29 and 37 for rotor and stator, respectively. The hub-tip ratio is 0.76. The design rotating speed is 17,000 rpm. Since this research is

focused on the subsonic conditions, the rotational speed is adopted as 12,000 rpm, so that the tip Mach number is less than 0.9. The design tip clearance is 98.7% of span. ,out

i,out in

=

out in

Compressor performance is obtained by the adjusting of the throttle area behind the stator. A flowmeter is equipped at the upstream of IGV to obtain flow rate. Total pressure is measured by two five-hole total pressure probes mounted at upstream and downstream of the compressor, respectively. Temperature probes are also used at both upstream and downstream to estimate isentropic efficiency. ,in / / 

The overall performance parameters are defined by the following equations. ,

$$\eta = \frac{H\_{\text{i,out}} - H\_{\text{in}}}{H\_{\text{out}} - H\_{\text{in}}} \tag{1}$$

$$
\pi = \frac{p\_{t, \text{out}}}{p\_{t, \text{in}}} \tag{2}
$$

$$\xi = \frac{\dot{m}\_{\text{in}}/\pi}{\dot{m}\_{\text{in}, \text{DP}}/\pi\_{\text{DP}}} \tag{3}$$

τ

where η is isentropic efficiency, *H*i,out isentropic outlet enthalpy, *H*in inlet enthalpy, *H*out outlet enthalpy, π total pressure ratio, *pt*,out isentropic outlet total pressure, *pt*,in inlet total pressure, ξ flow rate coefficient, . *<sup>m</sup>in* mass flow rate, . *min*,*DP* mass flow rate at design point, and π*DP* total pressure ratio at design point.

#### *2.2. Numerical Model*

The numerical simulation is accomplished by using a commercial CFD code, NUMECA. It is a popular solver for the flowfield simulation of turbomachinery. With the features of structured-grid supported, the finite-volume method and the suitable numerical treatments, the three-dimensional (3D) steady and unsteady flowfield can be obtained for the axial compressor.

The calculation area is between Section A and Section B in Figure 1. The upstream channel and downstream channel are extended as much as two IGV blade chord lengths and three stator blade chord lengths, respectively. Five cases with different tip clearance sizes from 0% to2% of chord length are established and marked as 0%, 0.5%, 1%, 1.5%, and 2% of the chord length of blade tip τ, respectively.

**Figure 1.** Sketch of the test rig.

Real air model is adopted as working fluid. The two-order central difference spatial scheme of Jameson is adopted in this study. Spalart–Allmaras (SA) model by Spalart et al. [24], which is widely applied in industry and academia, is also used in this study. For the inlet condition, total pressure, total temperature, and the velocity direction are specified, and refer to the experiment data. For the outlet condition, the static pressure is applied with the radial equilibrium. The rotating speed is applied on the rotor passage, while the other passages keep station. When processing the steady simulation, the mixing plane method is used to transfer the flowfield information between rotating

region and station region. Multigrid and local time stepping treatment are also applied to accelerate the calculation. The steady simulation is considered as converged when the mass flow is balanced and all of the residual lines are stable or fluctuate periodically. With the increasing of the outlet pressure, the last converged steady simulation is considered as near-stall (NS) condition. The unsteady method for this model had been discussed by Song et al. [25], including using the dual time step method of Jameson [26] and the sliding mesh method [27]. The 20 timesteps are set in every rotor passing period, and the number of the inner iteration step is set to 50. When the periodic fluctuation of the mass flow and the static pressure on the numerical probes is continued for ten passing times, it is considered as converged for the unsteady simulation.

#### *2.3. Grid Generation and Independent Test*

O4H-type topology is adopted in every row. The y+ is less than 5 to ensure capturing the blade loading more accurately. Mesh in tip clearance is refined to 17 layers in radial direction. Three sets of grids are generated to validate the grid independence. The grid number of every passage is shown in Table 1. From validation results shown in Figure 2, grid2 is appropriate for the calculation. A three level multi-grid is adopted to accelerate the calculation. The chosen grid is shown in Figure 3. In order to transfer the flow information between stational region and rotational region, the mixing plane method is used in this simulation.

**Table 1.** Grid number of three sets of grids (×10<sup>6</sup> ).


**Figure 2.** Validation for grid independence.

**Figure 3.** One passage of grid for simulation.

#### *2.4. Model Validation*

Based on the chosen grid, the performance curve is obtained by the numerical simulation. From Figure 4, by using the full-wheel simulation, the simulated mass flow range is enlarged and the numerical result is closer to the experiment data. It is found that the numerical results are always larger than experiment data; this might be the result of the neglection of the blade hub fillet in the numerical model. Even so, the numerical results match well with the experiment data on most mass flow conditions.

**Figure 4.** Performance comparison between experimental and numerical results.

#### **3. Results**

τ τ In this section, the influence of tip clearance on the performance is studied. The trend of performance change is discussed with the clearance increased from zero (0%τ) to large (2%τ), and the characteristics under peak efficiency (PE) working point and near stall (NS) working point are compared and discussed.

#### *3.1. Overall Performance*

Performance curves of every clearance model are shown in Figure 5. In a mass, the efficiency and pressure ratio are decreased as clearance size is increased, which is coherent with conclusions from predecessors' work. However, there is a certain tip clearance size, so that the stable working range can reaches to the maximum, which can be observed from stall margin change (Figure 6). The maximum surge margin is found in the design clearance model. When the tip clearance is larger or less than the design clearance model, the stall margin is decreased, obviously. In other words, with the tip leakage flow, the stability of the compressor might be enhanced. So, it is meaningful to find out the mechanism of leakage flow that turns the negative influence to positive influence on stall margin.

**Figure 5.** Performance curves of different tip clearancemodels. (**a**) pressure-flow rate; (**b**) efficiency-flow rate.

**Figure 6.** Stable working range under different clearance models.

The rotor loss factors along span is shown in Figure 7. The dimensionless rotor loss factor is defined in Equation (4).

$$F\_{\mathbf{r}} = \frac{P\_{i,t,r} - P\_{o,t,r}}{P\_{i,d,r} - P\_{o,d,r}} \tag{4}$$

where *P* is the pressure, subscript *i*, *o*, *t*, *d,* and *r* represent the rotor inlet, outlet, total, dynamic quantities and relative quantities, respectively. Along the span, the loss in tip region is the maximum. With the increasing of tip clearance, the loss in tip region is increased correspondingly. In addition, in tip region, loss factor is increased rapidly from 0%τ to 1%τ, which illustrates the negative impact of tip leakage on efficiency. When the tip clearance is over 1%τ, loss in tip region is stable with the increasing of tip clearance.

**Figure 7.** Loss factor distribution under different clearance models.

CP

Figure 8 shows the blade loading at 98% of span. Pressure coefficient is defined in Equation (5):

$$\text{CP} = \frac{P\_{\text{s}\rho} - P\_{\text{s},i}}{P\_{\text{s},i}} \tag{5}$$

where subscript *s* represents the static quantities. The blade loading of zero-clearance model is over two times larger than that of design-clearance. As the tip clearance is increased, the blade loading is decreased, which is benefit for the safety of the blade. This can be explained that the tip leakage flow can balance the pressure difference between the pressure side (PS) and the suction side (SS).

**Figure 8.** Pressure coefficient distribution at blade tip under different clearance models.

Ω *ω*

#### *3.2. Detail Results of Zero Clearance Model*

In this model, there is no tip leakage flow or tip clearance vortex at all. At 96% of span, relative Mach number contours with streamline of relative velocity are presented in Figure 9a. Under PE condition, there is seldom low-Mach region, except the boundary layer on the blade surface. From the streamline, since the incidence angle matches with the stagger angle, the flowfield is well without any vortex or separation. Under NS condition, there is a large area of the blockage region from the LE to the neighboring LE, and from the SS to PS. Since the mass flow is decreased, the incidence angle is far from the blade stagger angle, so a vortex is generated near the PS close to LE, and enrolls the surrounding fluid, with which the leading-edge spillage is occurred.

**Figure 9.** Time-averaged results at 96% of span. (**a**): Relative Mach contours with streamline of relative velocity. (**b**): Entropy contours (J/(kg K)).

Static entropy contour is often used to illustrate the loss distribution. Moreover, it is often used at tip span to locate the boundary of the tip leakage flow and the separation flow, because high loss will be induced when circumferential tip leakage flow/separation flow and main flow are interacted, in which the vortex would be generated and developed, and by which the entropy gradient is largest among tip-span contour. Static entropy contours are shown in Figure 9b. Under PE condition, the loss in the tip region consists of the wake loss from IGV, the boundary viscous loss and the wake from the rotor blade. The loss area and the amplitude are both low. Under the NS condition, it is found that the high-gradient interface is shifted forwardly to the front of the LE, which is coherent with the LE spillage from Figure 9a. Compared with the figure of PE, the loss region and the amplitude are

highly increased. The separation loss occupies the whole passage inlet, by which the loss on the whole passage, especially on the pressure, is highly increased.

In order to observe 3D distribution of separation and blockage, streamlines in the tip region are drawn in Figure 10 with velocity flux contours at passage cutting planes. The blockage area is attached to the casing-suction corner and reduced with the decreasing of radial distance. Along the suction side surface, streamline flows toward the tip region. When the flow arrives to casing, it turns to the middle passage along the casing side, and finally flows to the downstream due to the mixing process with the main flow. From the 3D streamline in Figure 10a, and limiting streamline in Figure 10b, the separation and blockage along the suction side is caused by leading-edge separation and the corner separation.

**Figure 10.** The three-dimensional (3D) figures for the passage flowfield. (**a**): streamlines in tip region with velocity flux contours at passage cutting planes. (**b**): streamlines in tip region with velocity flux contours at passage cutting planes.

For PE and NS conditions, here are the iso-surfaces of vorticity at υ = 0.25 in Figure 11. υ is defined as follows.

$$\upsilon = \frac{\Omega^2 + \mathcal{S}^2}{\omega^2} \tag{6}$$

where Ω, *S*, ω are the symmetric and antisymmetric parts of the gradient velocity tensor. Since the flow is very steady, only the steady result is shown for the PE condition.

**Figure 11.** Iso-surfaces of instantaneous vorticity at *λ*λ = 0.25. (**a**): Peak efficiency (PE), steady results. (**b**): Near-stall (NS), transient result, PBPF is the period of the blade passing.

It is different between the vortex distribution in PE and NS condition. For PE condition, there is no large vorticity cell in flow passage. However, there are two parts that are worthy to pay attention to. In the front and rear part of suction side, some little vorticity cells attach to the corner, from which the effect of corner separation is illustrated. In addition, vorticity cells in front of LE develop in a circumferential direction, which is the result of the effect of LE separation and casing separation. When the condition reaches to NS, from the black box in Figure 11b, these vorticity cells grow up to large cells, and shed with the period of blade passing. Flow incidence in front of the rotor passage will be increased because of the blockage, which will cause larger separation and corresponding larger blockage in the next passage. When the flow incidence increases to extent, leading edge spillage will occur, as shown in Figure 9a.

In Figure 12, the axial component of the velocity is shown in order to illustrate the unsteadiness of the velocity under NS condition. Although there is the vortex generation and moving by the leading-edge separation and corner flow, the unsteadiness is still very low. From the <sup>1</sup> 4 of the period of the blade passing (PBPF) to 4/4 PBPF, there is seldom unsteadiness that can be observed. This is because the blockage decreases the amplitude of the unsteadiness. In order to qualify the unsteadiness in NS condition, we set static pressure probes in front of LE, and the spectrum result is shown in Figure 13. The spectrum is dominated by character frequency related to rotation, including rotation frequency and blade passage frequency.

**Figure 12.** Axial component of velocity under NS condition, 99.35% of span, PBPF is the period of the blade passing.

**Figure 13.** Pressure spectrum in front of LE under zero clearance model.

#### *3.3. Detail Results of Design Clearance Model*

Figure 14 delineates distribution of the time-averaged axial momentum of tip leakage flow. MD represents the condition at which the flow rate is intermediate between PE and NS. µis defined as follows: *μ*

$$
\mu = \int\_{r\_{\rm tip}}^{r\_{\rm casing}} \frac{\rho V\_n V\_t}{\dot{m}\_{\rm in} V\_z / c\_x} \mathrm{d}r \tag{7}
$$

where . *min* is mass flow rate at rotor passage inlet, *c<sup>x</sup>* is chord length of blade, subscript *n* is normal direction, *t* is tangent direction, *z* is axial direction. 

**Figure 14.** Momentum distribution of tip clearance flow along chord.

μ From this distribution, the axial momentum is concentrated at the front part of the chord. As the flow rate decreases, the centroid of µ moves to the leading edge, hence, the leakage flow in the rear of the clearance is smaller. The effect of the existence and the distribution of leakage flow will be discussed in following paragraphs.

Figure 15 shows contours of the Mach number and entropy at 99.35% of span. In general, with the influence of the tip leakage flow, flow field in the tip region is much different with the zero-clearance model. From Mach contours, in both PE and NS conditions, flow separation and serious blockage are suppressed at the suction side, which is the result in mixing with circumferential leakage flow. However, in the NS condition, there is a block of a low speed region along the pressure side, which is due to the tip clearance vortex. The leading-edge spillage also occurred, but the reason is not the leading-edge separation, but the TCV. From entropy contours, the high-gradient interface has moved to the rotor inlet. Hence, the main flow could not flow pass the passage, which caused the low speed region. However, with the circumferential flow from the tip clearance, the low speed region will not develop to a large blockage. When it is at a rear part of the passage, the low speed region could cross the passage and expend to the suction side. This is because the tip leakage flow in the rear part of the clearance is too small to supply enough momentum for a low speed region. As the flow in the tip region is dominated by the tip clearance vortex and tip leakage flow, from the streamlines at the passage outlet in Figure 16, in the rear part of clearance, there is the vortex generated by the sheared tip leakage flow.

**Figure 15.** Time-averaged results at 96% of span. (**a**): Relative Mach contours with streamline of relative velocity. (**b**): Entropy contours (J/(kg K)).

**Figure 16.** Streamline on suction side, casing, and passage outlet.

Figure 17 shows the axial component of the velocity, which can represent the location of the clearance vortex. PTCV is a period of the shedding of the tip clearance vortex, which corresponds to about 0.5 Blade passing frequency (BPF). The tip clearance vortex results in the leading-edge spillage and mixing process between the main flow and tip leakage flow. Tip clearance vortex is formed near the leading edge, developed along the rotor passage inlet, interacted with the leading edge of the other side, and reduced along the pressure side. Streamlines from the tip clearance are shown in Figure 18, which can show the detail of the tip clearance vortex. Streamlines from the tip clearance roll up to a vortex, which is the result in the mixing process, with main flow near the leading edge. As the axial momentum ratio between the leakage flow and main flow increases, the vortex trajectory moves to the rotor inlet. When the streamlines reach the next leading edge, they are diverged. Some streamlines cross in front of the leading edge to the next passage, which is called leading-edge spillage. Streamline in the core of the vortex is developed and finally stagnated. Then the vortex breakdown occurs (as yellow blank locates), because the flow is a counter-pressure along the vortex trajectory. It can trigger several low speed cells (as the red dot line shows); some of them will be developed to become larger at the mid-passage. The streamlines near the large low-velocity cell begin to diverge by the blockage effect of low-velocity cells, which is called secondary vortex as black arrow shows. Moreover, when the larger low-velocity cell is developed near the next leading edge, the spillage streamlines will be increased by the blockage effect.

**Figure 17.** Unsteady contours of axial component of velocity. (PTCV is the period of the TCV (tip clearance vortex) shedding).

**Figure 18.** Unsteady streamlines from tip clearance with the iso-surface of Ma<sup>r</sup> = 0.2, PTCV is the period of the TCV (tip clearance vortex) shedding.

In addition, the unsteadiness in the tip region is enhanced (a lot) by vortex breakdown and secondary vortex, which is illustrated in Figure 19. Compared with the spectrum in PE and NS condition, the amplitude increasement at BPFs are limited. However, the amplitudes at the frequencies related to the tip clearance vortex and rotation are increased from nearly zero. Compared with the spectrum of zero-clearance model in the NS condition, amplitude at BPFs are increased, while amplitudes at rotation frequencies are about the same.

**Figure 19.** Pressure spectrum in front of LE under design clearance model. (**a**): under PE condition. (**b**): under NS condition.

To sum up, under design clearance model, the influence of tip clearance dominates flow field in the tip region. The separation at leading edge and casing-suction can be eliminated by the tip leakage flow. However, the low-velocity region can be triggered by the tip clearance vortex in NS condition. Moreover, the unsteadiness in the tip region is enhanced by the evolution of the tip vortex.

#### *3.4. Detailed Results of Large Clearance Model*

As the clearance is increased, the total leakage flow is increased correspondingly. From the Mach number in Figure 20a, the separation in the tip region is inhibited furtherly. However, as the tip leakage flow is increased in scale and value, flow structure in the tip region will be influenced, and the loss is increased with the tip clearance in Figure 20b. On the one hand, the radial scale of leakage flow is increased, which will enlarge the scale of the tip clearance vortex. On the other hand, the axial momentum ratio between leakage flow and main flow is increased, which is the reason for the forward-moving of the tip clearance vortex. Therefore, compared with design clearance, the forward spillage occurred in a larger flow rate. In addition, the low-velocity cell can only be developed to the front–middle passage, then its velocity will increase by the increased circumferential leakage flow. The unsteadiness in the tip region is also enhanced by the increasing of the tip clearance, which is illustrated by the moving of the low-velocity cell in Figure 21. As the spectrum shows in Figure 22, compared with that from the zero-clearance and design-clearance models, the amplitude at BPFs increased a lot.

**Figure 20.** Time-averaged results at 96% of span. (**a**): Relative Mach contours with streamline of relative velocity. (**b**): Entropy contours (J/(kg K)).

**Figure 21.** Unsteady streamlines from tip clearance with the iso-surface of Ma<sup>r</sup> = 0.2, PTCV is the period of the TCV (tip clearance vortex) shedding.

**Figure 22.** Pressure spectrum in front of LE under large clearance model.

#### **4. Conclusions**

This paper studies the influence of tip clearance on flow characteristics related to performance. Based on full-passage grids, the validated results are obtained by the means of the RANS method. A series of performance curves are obtained in different clearance models, from which it is found there exists an optimum clearance for a stable working range. It is found that the tip leakage flow can benefit the stable working range if the tip clearance is small enough, while it turns out to be harmful to the stability if the tip clearance is too large.

Under the zero clearance model, it is the separation at the casing-suction corner that dominates the flow field in the tip region. Under this blockage, incidence of upstream is increased, which will trigger larger separation in the next passage. Finally, some passages are fully blocked and stall occurs. As the tip clearance is increased from zero to the optimum value, the separation and blockage along the suction side are suppressed by tip leakage flow.

Under the design clearance model, the influence of the tip clearance dominates flow field in the tip region. The existence of the tip leakage flow can benefit to the reduction of separation at the leading edge and casing-suction. However, the low-velocity region can be triggered by the tip clearance vortex in NS condition. Moreover, the unsteadiness in the tip region is enhanced by the evolution of the tip vortex.

When the tip clearance increases to extent the tip leakage flow increases in scale and value. Therefore, the scale and unsteadiness of the tip clearance vortex increases correspondingly, which can trigger the forward spillage earlier.

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

**Funding:** This research was funded by the National Science and Technology Major Project (2017-II-0006-0019).

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

#### **Nomenclature**


#### **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* **The Influence of the Blade Outlet Angle on the Flow Field and Pressure Pulsation in a Centrifugal Fan**

#### **Hongchang Ding \*, Tao Chang and Fanyun Lin**

College of Mechanical and Electronic Engineering, Shandong University of Science and Technology, Qingdao 266590, China; 17863935466@163.com (T.C.); 15662728726@163.com (F.L.)

**\*** Correspondence: dhchang@sdust.edu.cn; Tel.: +86-178-5272-8669

Received: 29 September 2020; Accepted: 5 November 2020; Published: 8 November 2020 -

**Abstract:** This paper takes centrifugal fan as the research object and establishes five impeller models with different blade outlet angles. By means of computational fluid dynamics (CFD), the external characteristics of the centrifugal fan and the internal characteristics, including the velocity, pressure, and turbulent energy distribution, at the middle span plane of the impeller or fan were obtained and compared. In addition, the pressure fluctuations surrounding the impeller outlet were also analyzed. The results showed that the change of the blade outlet angle of the centrifugal fan had a great influence on the performance; the total pressure and efficiency of the fan were the highest when the outlet angle of the blade was increased to 29.5◦ under the design flow rate; and the influence of the outlet angle on the fan performance was different in off-design conditions. On the other hand, at different flow rates, the change of the internal flow field with the increase of the outlet angle was different. For the pressure fluctuation of the fan, by increasing the blade outlet angle properly under high flow conditions, the fluctuation amplitude of the fan at the blade frequency and its frequency multiplication could be reduced, which is conducive to decreasing the impeller noise. The research results have good guiding significance regarding the design of the pneumatic performance and noise reduction performance of centrifugal fans.

**Keywords:** centrifugal fan; blade outlet angle; aerodynamic performance; numerical simulation

#### **1. Introduction**

Fans belong to the category of general machinery and are widely used in various industries of national economy. They are indispensable equipment for industrial and agricultural production. According to statistics, the power consumption of wind turbines accounts for 8–10% of the total power generation in China (Chen [1]). At present, centrifugal fans occupy a large proportion in China's energy system. Therefore, it is of great significance to research and improve the pneumatic performance of fans for energy saving. However, in the process of optimizing the performance of fans, the traditional experimental methods have long cycles and high cost, and it is difficult to visually display the gas distribution inside a fan. Therefore, a CFD (Computational Fluid Dynamics) numerical simulation technology that can effectively reduce the sample size of experimental design and capture the details of flow inside the fan more specifically, has increasingly caught the interest of scholars. Many scholars have used numerical simulation technology to study centrifugal fans (Zhang et al. [2]; Zhou et al. [3]; Lin et al. [4]; Yu et al. [5]; Kishokanna et al. [6]).

The impeller is the main moving part of a centrifugal fan, and the structural parameters of the impeller include the blade shape, blade profile, outlet width, number of blades, inlet and outlet diameter, etc. An excellent impeller design is helpful to improve the aerodynamic performance of the fan. At present, many scholars have studied the influence of certain impeller parameters on the aerodynamic performance of fans.

Li et al. [7] studied the influence of blade shapes on the performance of high specific speed centrifugal fans, and found that the blockage phenomenon at the blade outlet of impeller with plate blade and the turbulent kinetic energy inside the volute were weakened under the condition of a large flow rate, so that the performance of the fan with plate blade was better than that of the airfoil blade.

Wu et al. [8] compared the performance of a centrifugal fan with different blade profiles, and found that the centrifugal fan with a double-arc blade was higher in the efficiency and total pressure under the design condition, but the axial power consumed by equal deceleration blade was smaller. However, under the condition that other parameters of the fan remain unchanged, the internal flow of the fan with an equal deceleration blade was more uniform under the condition of low flow.

Jian et al. [9] found that when the blade outlet width changed, various losses of the fan increased, and the efficiency decreased. With the decrease of the blade outlet width, the flow-pressure curve of the fan shifted to the lower left and the pressure decreased with the increase of the flow rate. This provides a reference for the design of the outlet width of the impeller and the reconstruction of the impeller. Liu et al. [10] found that the aerodynamic performance of the fan can be improved by increasing the number of blades and the diameter of the blade outlet. The optimized fan with a 12-blade number and increased blade outlet diameter was better than the prototype fan in terms of the total pressure and efficiency.

In addition, some scholars studied the interaction among several structural parameters of the impeller. For example, Esra et al. [11] used a neural network method to determine the optimal family of the impeller inlet and outlet radius and inlet and outlet angle to reduce the noise level in the early stage of fan design. Meng et al. [12] studied the influence of three blade structure parameters on the performance of a centrifugal fan, which were inlet angle, outlet angle, and blade number, and obtained that the best combination of the three blade structure parameters using the response surface model (RSM) optimization method. The maximum efficiency was 93.7%.

Shi et al. [13] optimized the combination of the impeller and guide vane by studying the internal flow law of fans and combining this with numerical calculation results. They found that the streamline in the optimized guide vane was uniform and that the vorticity was reduced compared with the original guide vane. HEO et al. [14] analyzed the aerodynamic characteristics of a centrifugal fan with additional splitter blades in the impeller by using three-dimensional Reynolds-averaged Navier–Stokes (RANS). The global Pareto optimal frontier for centrifugal fan design was obtained by using a hybrid multi-objective evolutionary algorithm and response surface approximation model.

As the main parameter of impeller, the blade outlet angle of centrifugal fans has also been studied by scholars. Wang et al. [15] optimized the fan design based on an orthogonal experimental design method and CFD numerical simulation and obtained the relative optimal combination model of the impeller geometric parameters and speed. Through the range analysis of the calculation results, they found that the blade outlet angle had the greatest impact on the fan efficiency.

Recently, Swe et al. [16] discussed the flow characteristics of centrifugal fans with different blade outlet angles using the CFD numerical simulation method. The results showed that when the blade outlet angle was 25◦ , the steady flow and unsteady flow were more uniform, and the flow distribution in the circumferential direction had little change. Yu et al. [17] also studied the effect of the blade outlet angle on the performance of a multi-blade centrifugal fan, and found that the wind pressure and efficiency of the fan increased in the flow range of 420–725 m<sup>3</sup> /h with the increase of the blade exit angle, and properly increasing the blade outlet angle can reduce the pulsation amplitude of the fan at the blade frequency and its multiplier frequency.

Although many scholars have studied the blade outlet angle on the centrifugal fan, few have studied the effect of the blade outlet angle on the internal flow field and pressure fluctuation of the straight blade fan. In this paper, a type of straight plate blade fan was selected as the research model. Then, on the basis of the original blade outlet angle, five outlet angles (26.5◦ , 28◦ , 29.5◦ , 31◦ , and 32.5◦ ) were obtained, and the geometric model was established using SolidWorks software. ANSYS CFX software was used for the numerical simulation, and the feasibility of the calculation results was

verified by experiments. Through the software post-processing, the influence of different blade outlet angles on the internal flow field and performance of the centrifugal fan was studied. The amplitudes of fan blades with different outlet angles were obtained and analyzed using fast Fourier transform. The research results can provide certain reference values for the efficient, safe, and stable operation of a centrifugal fan.

#### **2. Research Model and Simulation Method**

#### *2.1. Research Model*

This paper took a certain type of centrifugal fan as the research model, and the main design parameters of the fan are shown in Table 1.


**Table 1.** Main design parameters of the centrifugal fan.

The strategy of this paper was to change the blade outlet angle of the original fan model, and its specific impeller parameters are shown in Table 2. The overall geometric model was established in SolidWorks as shown in Figure 1.

**Figure 1.** Geometric model of the original fan.



#### *2.2. Mesh and Check of Grid Independence*

δ In this paper, SolidWorks was used for the three-dimensional modeling of the centrifugal fan. In the process of modeling, the geometric model of the original fan was simplified, and some unimportant chamfers, fillets, and gaps were ignored. Based on this, the flow-path model of the centrifugal fan was established. The flow-path model was divided into three parts: air inlet flow-path, impeller flow-path

β

and volute flow-path, and they were combined to form the centrifugal fan model. To ensure the real working condition of the air flow, the inlet and outlet flow-path model were extended properly.

After completing the establishment of the fan model, the mesh generation of the model began. Due to the simplicity of the axisymmetric model of the fan inlet section, the hexahedron structure was used for mesh generation. However, the impeller and volute were meshed with the tetrahedron structure due to the complexity of their flow passages, and the meshing of impeller and volute is shown in Figure 2. To make the numerical simulation as close to the actual situation as possible, all boundary layer grids were refined to meet the requirements of the wall function method. The minimum orthogonal angle of the mesh was greater than 27◦ , and the maximum extension ratio was less than 2.1, which will ensure high mesh quality.

**Figure 2.** Mesh generation of the impeller and volute. (**a**) Impeller. (**b**) Volute. (**c**) Local refinements.

Table 3 shows the grid independence verification results of the centrifugal fan at design points, in which the total pressure and efficiency were obtained under five different numbers of grids. When the number of grids increased from 3,284,561 to 495,126, the total pressure and efficiency tended to be stable. Based on the accuracy of the simulation results and the cost of the computing time, the scheme of No. III was finally adopted.


**Table 3.** Scheme for grid independence check.

#### *2.3. Simulation Settings*

ε At present, the standard k-ε turbulence model is the most widely used turbulence model, which has good stability and fast calculation speed compared with the zero equation model and the single equation model, and it is very suitable as the research object of this paper.

The numerical simulation was carried out in ANSYS CFX14.5, and the simulation process was based on the following settings:


To study the possible influence of different blade outlet angles on the pressure fluctuation characteristics of the impeller outlet and volute of fan, monitoring points were set on the outlet of the impeller with mid-section (z = 79 mm) to monitor the unsteady flow. The layout of the monitoring points is shown in Figure 3, and the impellers with different blade outlet angles had the same monitoring position. In Figure 3b, P1 to P4 are distributed clockwise along the circumference, and P5 and P6 are arranged in the volute tongue and volute outlet section, respectively.

**Figure 3.** Schematic: (**a**) mid-section (z = 79 mm); (**b**) location of the monitoring point.

#### **3. Experimental Verification**

#### *3.1. Experimental Equipment*

To verify the feasibility of the established three-dimensional model of a centrifugal fan in numerical simulation, the external characteristics experiment of the fan model was carried out. The fan test was conducted according to the GB / T1236–2000 standard, the air intake test system was adopted, and the test device was type C (pipe inlet and free outlet) specified by the national standard. The schematic diagram of the experimental device is shown in Figure 4.

**Figure 4.** Schematic diagram of centrifugal fan experimental device.

The centrifugal fan test platform was primarily composed of the fan inlet pipe, main body of the centrifugal fan, and motor. Among them, the inlet pipe part mainly included the current collector, compensation micromanometer, rectifier network, and U-tube liquid pressure gauge. Table 4 shows the instrument used in the test and its specific parameters. The actual drawing of the fan is shown in Figure 5, and Figure 1 shows the geometric model of Figure 5.



≤

≤

°

**Figure 5.** The actual drawing of the fan.

#### *3.2. Experimental Results*

Figure 6 shows the results comparison between the numerical simulation and experiment. The longitudinal coordinate is the total pressure coefficient *P* and total pressure efficiency η, and the calculation equation is expressed as:

$$
\overline{P} = \frac{P}{\rho u^2},
\tag{1}
$$

$$
\mu = \frac{\pi D\_2 n}{60},
\tag{2}
$$

$$
\eta = \frac{Q P}{M \omega'} \tag{3}
$$

where *P*, ρ, *u*, *n*, *Q*, *M* and ω respectively represent the total pressure, Pa; air density, kg/m<sup>3</sup> ; impeller outlet diameter, mm; rotational speed, r/min; experimental flow rate m<sup>3</sup> /s; torque, N·m; angular velocity, rad/s.

From Figure 6, the experiment results are close to the numerical results. Q/Qdes is the ratio of the experimental flow rate to the design flow rate. As the numerical simulation simplifies the model and neglects the loss, the simulation value was higher than the experimental value. The error between the tested and calculated total pressure was the smallest near the work point but gradually increased with the flow away from the design point, and the maximum deviation value did not exceed 8%.

There were two main influencing factors: One was the complex internal flow rule of fans under a small flow rate; the turbulence model selected in this paper could not effectively simulate the flow separation in the flow channel due to time and condition limitations. The other factor was that the actual flow rate of the fan was transient, which is ignored by the frozen rotor method. Although there were some errors between the experimental and the numerical results, the general changing trend of the total pressure and efficiency curves was the same, and the reliability of the numerical simulation can be confirmed according to the experimental results.

**Figure 6.** Curve of the total pressure coefficient and efficiency of the centrifugal fan model.

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

#### *4.1. Contrast of External Characteristics*

β From Figure 7, the total pressure of the fan decreased with the increase of the flow rate. The efficiency of the fan increased first and then decreased, and it reached the highest value near the rated operating point. At Q/Qdes < 1.4, the total pressure of the fan increased first and then decreased with the rise of the blade outlet angle at the same flow rate. When β<sup>2</sup> = 29.5◦ , the total pressure was the highest. In addition, the change trend of efficiency was different from the total pressure change. The outlet angle of the blade with the highest efficiency was different under different flow rates. At a

high flow rate, the efficiency increased with the increase of the outlet angle. At the design or a low flow rate, the efficiency was the highest when the blade outlet angle was 26.5◦ and reached the maximum of 84.85% at the designed flow rate.

β

**Figure 7.** External characteristic curve of fan with different blade outlet angles: (**a**) Comparison of the total pressure coefficient. (**b**) Efficiency comparison.

#### *4.2. Velocity Vector Distribution*

To further describe the internal flow of the fan, the relative velocity distributions in the middle section of the flow passage were obtained, as shown in Figure 8. There were vortices of different degrees in the impeller channel between 90◦ and 270◦ (see Figure 3) at different flow rates, and the number of vortices gradually decreased with the increase of the flow rate. At a low flow rate, the outlet angle had little effect on the relative velocity distribution; however, the velocity near the volute tongue increased significantly with the increase of the outlet angle. β

**Figure 8.** Velocity vector distributions at the middle span plane of the fan.

At the design flow rate, when the blade outlet angle was 31◦ , the number of impeller channels with a vortex was the least, and the low-speed area of the impeller outlet was the least. Under the condition of a large flow rate, there were low-speed regions of varying degrees near the impeller outlet between 90◦ and 270◦ , and the area of the low-speed area decreased with the increase of the blade outlet angle. However, the flow state of the impeller channel near the volute tongue became increasingly unstable. When β<sup>2</sup> = 32.5◦ , the low velocity region and flow separation appeared on the suction surface of the trailing edge of the blade, resulting in flow blockage. It may be that the outlet angle was too large to match the volute and that the fluid impact on the volute tongue caused this phenomenon.

#### *4.3. Pressure Distribution*

−

Figure 9 shows the static pressure distribution of the middle span plane of impellers with different blade outlet angles and different flow rates. Under the condition of a low flow rate, the low-pressure region of the impeller inlet increased first and then decreased with the increase of the outlet angle. At the design flow rate, the area of the low-pressure zone increased first and then decreased with the increase of the impeller outlet angle between 90◦ and 270◦ . Under a large flow rate, there was a reverse pressure gradient in the impeller channel between 90◦ and 270◦ , and this phenomenon was alleviated with the increase of the blade outlet angle.

The static pressure increased with the increase of the flow rate. The maximum static pressure occurred at the impeller outlet near the volute tongue, which increased with the increase of the blade outlet angle, and the area of high pressure area became larger. In the static pressure distribution diagram of the impeller, we observed that the static pressure value of the fan gradually increased due to the continuous work of the impeller blade rotation on the air flow. At the impeller inlet, there was a clear low pressure area on the suction surface of the blade, and reverse flow occurred near the wall, resulting in separation loss. At high flow rates, this phenomenon was alleviated with the increase of the outlet angle.

**Figure 9.** Static pressure distributions at the middle span plane of the impeller.

As can be seen from Figure 10, the total pressure distribution of the fan gradually increased from the impeller inlet to the outlet. As the flow rate increased, the total pressure increased. Under different outlet angles, the pressure changes in the impeller passage mainly occurred between 90◦ and 270◦ . At a low flow rate, the total pressure did not change greatly with the increase of the blade outlet angle. At the design flow rate, the low-pressure region appeared at the outlet of the single impeller, which increased first and then decreased with the increase of the outlet angle.

Combined with the velocity vector diagram, the low-pressure area appeared in the low speed region. At high flow rates, the total pressure increased with the increasing of the outlet angle. With the increase of the blade outlet angle, the effective working area of the impeller increased, the total impeller pressure increased overall, especially in the circumferential area of the outlet of the blade. In addition, the pressure distribution of different flow paths between blades was different in the range of 90◦−270◦ , which indicates that the flow characteristics inside the fan were asymmetric.

**Figure 10.** Total pressure distributions at the middle span plane of the impeller.

#### *4.4. Turbulence Kinetic Energy Distribution*

It can be seen from Figure 11 that the high turbulent kinetic energy region was mainly distributed near the volute tongue at the impeller outlet and gradually diffused toward the surrounding area. When the gas flowed through this area, boundary layer separation easily occurred. In general, the turbulent energy was the smallest when the blade outlet angle was 29.5◦ . With the increase of the blade outlet angle, the turbulent energy of the volute section first decreased and then increased. The outlet angle of the blade had a suitable value, which can reduce the flow loss of the volute section and improve the efficiency of the fan.

**Figure 11.** Turbulence kinetic energy distributions at the middle span plane of a fan.

#### *4.5. Pressure Pulsation Analysis*

To better analyze the influence of the blade outlet angle on the pressure pulsation in the impeller outlet area, the time-domain characteristics of pressure pulsation at each monitoring point of the impeller were transformed by fast Fourier transform, and the frequency domain of pressure pulsation was obtained, as shown in Figure 12. It can be seen from the diagram that the main frequency of pulsation at the monitoring point of the impeller outlet of the centrifugal fan was the blade passing frequency (fBPF) and its frequency multiplication, and the amplitude of pulsation reached the maximum at the blade passing frequency. Due to rotor–stator interactions, pressure fluctuations near the volute tongue (P1, P5) were large. However, from the view of the frequency domain map of the monitoring points, the change rules of each monitoring point were different at different outlet angles, so further analysis is required. The blade passing frequency was calculated by Equation (4).

$$f\_{\rm BPF} = \frac{\rm NZ}{60} \,\prime \tag{4}$$

where N is the speed, r/min; Z is the number of blades.

For further analysis, we use mean pressure amplitude *c<sup>p</sup>* to calculate pressure pulsation energy of all the six measuring points for different outlet angles.

$$
\overline{c\_p} = \frac{\sum\_{1}^{n} c\_{p-i}}{\mathsf{6}} (n = \mathsf{6}),
\tag{5}
$$

where *<sup>c</sup>p*−*<sup>i</sup>* represents the pressure amplitudes at *<sup>f</sup>BEF* for different measuring positions.

− −

( ) **<sup>1</sup>**

**6**

**6**

**Figure 12.** Frequency domain of the pressure pulsation for different blade outlet angles: (**a**) 26.5◦ ; (**b**) 28◦ ; (**c**) 29.5◦ ; (**d**) 31◦ ; and (**e**) 32.5◦ .

Figure 13 shows the average pressure amplitude at different outlet angles. In general, with the increase of the blade outlet angle, the pressure fluctuation amplitude of the fan increased. However, the amplitude of pressure fluctuation decreased with the increase of the blade outlet angle at 26.5◦−28◦ and 29.5◦−31◦ . Therefore, overall, the outlet angle of the blade increasing in a certain range reduced the pressure pulsation on the impeller, which is beneficial to reduce the noise of the impeller.

**Figure 13.** Average pressure amplitude at different angles.

#### **5. Conclusions**

In this paper, five models of centrifugal fan impellers with different blade outlet angles were established, and the influence of the blade outlet angle on centrifugal fan performance was studied using CFD software. The conclusions can be drawn as follows:


**Author Contributions:** Conceptualization, H.D.; Data curation, T.C.; Methodology, H.D. and T.C.; Software, F.L.; Supervision, H.D.; Validation, F.L.; Writing–original draft, T.C.; Writing–review and editing, H.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been supported by key research and development project of Shandong province (2017GGX203005), Natural Science Foundation Project of Shandong province (ZR2019MEE068) and Design and performance optimization of magnetic levitation high speed centrifugal pump(201810424020). The supports are gratefully acknowledged.

**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* **A CFD-Based Shape Design Optimization Process of Fixed Flow Passages in a Francis Hydro Turbine**

### **Ujjwal Shrestha <sup>1</sup> and Young-Do Choi 2,\***


Received: 4 October 2020; Accepted: 27 October 2020; Published: 31 October 2020

**Abstract:** In recent times, optimization began to be popular in the turbomachinery field. The development of computational fluid dynamics (CFD) analysis and optimization technology provides the opportunity to maximize the performance of hydro turbines. The optimization techniques are focused mainly on the rotating components (runner and guide vane) of the hydro turbines. Meanwhile, fixed flow passages (stay vane, casing, and draft tube) are essential parts for the proper flow uniformity in the hydro turbines. The suppression of flow instabilities in the fixed flow passages is an inevitable process to ensure the power plant safety by the reduction of vortex-induced vibration and pressure pulsation in the hydro turbines. In this study, a CFD-based shape design optimization process is proposed with response surface methodology (RSM) to improve the flow uniformity in the fixed flow passages of a Francis hydro turbine model. The internal flow behaviors were compared between the initial and optimal shapes of the stay vane, casing, and the draft tube with J-Groove. The optimal shape design process for the fixed flow passages proved its remarkable effects on the improvement of flow uniformity in the Francis hydro turbine.

**Keywords:** CFD; shape optimization; Francis turbine; fixed flow passage; flow uniformity

#### **1. Introduction**

Hydropower is considered a reliable renewable source for electricity production. The hydraulic turbine is an essential component of the hydropower plant. Among various types of hydraulic turbines, Francis turbines are widely used over a wide range of flow rates and heads [1]. The main hydro passage parts of the Francis turbine are composed of a spiral casing, stay vane, guide vane, runner and draft tube. The flow instabilities in the fixed flow passages can cause failure in the whole hydro turbine system. The design of the fixed flow passages is dependent on the moving components of the runner and guide vane in the Francis turbine.

The fixed flow passages are designed to keep the proper flow uniformity by suppressing the pressure pulsation, vortex-induced vibration and swirl flow. The main objective of the stay vane is to maintain the uniform flow from the casing to guide vane and runner flow passages [2]. The non-uniform flow distribution from the stay vane causes the vortex-induced vibration, which initiates the failure in the stay vane [3]. The purpose of the spiral casing is to direct the fluid from the penstock pipe to the stay vane and guide vane. Kurokawa and Nagahara [4] explained the free-vortex, accelerating and decelerating types of the spiral casing. The flow behavior is dependent on the casing shape. The improper flow distribution causes pressure pulsation and secondary vortex, which induces the cracking in the casing wall. Price indicated that the severe pressure fluctuation in the spiral casing causes the brittle crack in the casing wall [5]. The draft tube is designed to improve the dynamic

energy in the runner outlet and recover the suction head [6]. The existence of swirl flow in the draft tube causes the flow instabilities [7–9]. J-Groove installation suppresses the flow instabilities in the draft tube of the Francis hydro turbine [10]. J-Groove is the groove engraved on the wall of the draft tube that induces reverse jet flow through the shallow groove channels to suppress the swirl flow [11].

The computational fluid dynamics (CFD) analysis has become one of the main tools for turbomachinery flow analysis. The application of CFD analysis makes it easier to evaluate a large number of design cases with precise and accurate results. It is used to predict the internal flow behavior of the turbomachinery, flow separation, and loss distribution in Francis turbine components. Many researchers have conducted CFD analysis on the Francis hydro turbine for the prediction of performance [12], part-load performance [13], suction performance [14], unsteady flow behavior [15]. The CFD analysis and optimization techniques were integrated for the optimization of the runner blade [1]. The maximum improvement in the moving components (runner and guide vane) was achieved via CFD-based optimization [16–19]. However, few studies related to the fixed flow passages (spiral casing, stay vane and draft tube) of Francis hydro turbines are available [20,21].

Nowadays, design optimization using numerical analysis is widely used for turbomachinery. Wu et al. performed the CFD-based design optimization for a Francis hydro turbine. They showed a comparison between the initial and optimal design of turbines at the design point [19]. They mainly focused on the optimization of the runner blade of the Francis hydro turbine. The conventional blade design approaches integrated with the advanced CFD analysis are powerful and effective tools for the design optimization of turbomachinery. A CFD-based design optimization system that integrates internally developed parametrized mathematical geometry models, automatic mesh generators and commercial 3D Navier–Stokes code like ANSYS CFX 19.2 permits designers to interactively generate, modify and visualize the geometric model of turbine components. The design process can be repeated until a fully optimized model with satisfactory performance is obtained. Nakamura and Kurosawa [22] conducted the design optimization of a high specific speed Francis turbine using a multi-objective genetic algorithm (MOGA). The design optimization of hydraulic machinery can be performed by multilevel CFD techniques [23]. The multilevel CFD technique reduced the computation time. Sosa et al. [6] performed the design optimization of the draft tube by using CFD analysis. Si et al. [24] proposed a multi-point design process based on CFD and an intelligent optimization method for the automotive electronic pump. Ayancik et al. [1] conducted a simulation-based design of a Francis hydro turbine runner that was performed by following a surrogate model-based optimization. The conventional CFD-based design process is executed through trial and error; hence, designing a runner for a Francis hydro turbine can take several months. Due to these drawbacks of conventional CFD-based design, CFD-based optimization design approaches can be followed for the reduction of calculation time and better shape design.

It is essential to integrate a robust and flexible design tool in a CFD-based design optimization system to allow automatic generation and modification of the design geometry. The objective of this study is to propose a CFD-based shape design optimization process for fixed flow passages in the Francis hydro turbine. For the CFD-based shape design optimization, the surrogate model was prepared by using response surface methodology (RSM). The various RSMs were evaluated by the goodness of fit test for the precise and accurate response surface. The multi-objective genetic algorithm (MOGA) was applied for the optimization of the fixed flow passages. The hydraulic design and optimization framework for stay vane, casing, and a draft tube with J-Groove can be generalized for reaction hydro turbines (Francis turbine and Pump turbine). The CFD-based optimization process included the parametric design of stay vane, casing, and the draft tube with J-Groove, fluid domain modeling, meshing, ANSYS CFX solver, post-processing, design of experiment (DOE), response surface preparation and multi-objective optimization.

#### **2. Hydro Turbine Design and Optimization Methodology**

#### *2.1. Hydro Turbine Design*

#### 2.1.1. Process of Hydraulic Design

Figure 1 illustrates the hydraulic design process of the Francis hydro turbine proposed in this study. The conceptual design of the turbine components is based on the turbomachinery theory. The conceptual design for the Francis hydro turbine was prepared according to the turbine specification. Generally, the hydro turbine design is commenced with the runner design. The guide vane is designed according to the flow angle at the runner inlet. The stay vane design should match the inlet flow angle of the guide vane. The proper flow distribution at the stay vane inlet should be maintained by casing design. The runner outlet flow angle is a constraint for the draft tube design. A fixed flow passage design is linked with each other. Therefore, the initial shape design was completed in serial order as in the conceptual design by theory and detail design by 3D shape modeling. It was challenging to obtain the whole turbine passages optimization at once because it consisted of numerous design variables and overlapping constraints. Hence, the fixed flow passages shape was optimized for each passage separately to reduce the computational cost and make an effective optimization process. In this study, runner and guide vane design conditions were fixed, which created the constraints for the optimization process followed.

**Figure 1.** Hydraulic design process of the Francis hydro turbine with computational fluid dynamics (CFD)-based optimization.

2.1.2. Francis Turbine Specification and Performance

݊√ܲ The design specification of the 100 MW class Francis hydro turbine model is shown in Table 1. The minimum and maximum heads of the Francis hydro turbine are 66.5 m and 110 m, respectively. The design flow rate of the prototype turbine is 125.4 m<sup>3</sup> /s. The turbine maximum output power is 113 MW, and the minimum output power is 62.3 MW. The specific speed *N<sup>s</sup>* , unit discharge *Q*<sup>11</sup> and unit speed *N*<sup>11</sup> are evaluated by using Equations (1)–(3), respectively.

ܰ<sup>௦</sup> =

ହ

$$N\_s = \frac{n\sqrt{P}}{H\_\sqrt^5} \tag{1}$$

$$Q\_{11} = \frac{Q}{D\_\varepsilon^2 \sqrt{H}}\tag{2}$$

$$N\_{11} = \frac{nD\_\varepsilon}{\sqrt{H}}\tag{3}$$

where *n* is the rotational speed in min−<sup>1</sup> , *P* is the output power in kW, *H* is effective head in m, *Q* is flow rate in m<sup>3</sup> /s, *D<sup>e</sup>* is the runner outlet diameter.


**Table 1.** Design specification of the 100 MW class Francis turbine.

Francis turbine fluid domain by initial design is shown in Figure 2a, which consists of a spiral casing, 20 stay vanes, and 20 guide vanes with 13 runner blades, and the elbow-type draft tube. The runner inlet and outlet diameters are *D<sup>i</sup>* = 4863 mm and *D<sup>e</sup>* = 3995 mm, respectively. The installed capacity of the Francis turbine is 100 MW. The efficiency hill chart of the Francis hydro turbine by initial design and CFD analysis is shown in Figure 2b. The various guide vane openings are used to regulate the flow rate for the Francis hydro turbine. The guide vane opening from 8◦ to 41◦ is used to change the flow rate. The design point for the Francis hydro turbine is determined at *N*<sup>11</sup> = 76.12 and *Q*<sup>11</sup> = 0.87. The best efficiency of the Francis hydro turbine is located in the range of guide vane angles of 23◦ to 26◦ and unit speeds of *N*<sup>11</sup> = 70 to *N*<sup>11</sup> = 80. − −

**Figure 2.** (**a**) Schematic view and (**b**) efficiency hill chart of the 100 MW class Francis hydro turbine by initial design.

#### 2.1.3. Stay Vane Design

*δ*

The purpose of the stay vane (SV) is to guide the water flow from the casing to guide vane and runner, and for structural purpose. The main design parameters for the stay vane are vane angles,

> , ܽ, ܽ௧ሿ ்

<sup>ହ</sup>ߙ,..., <sup>ଵ</sup>ߜ , <sup>ହ</sup>ߜ,...,

ௌ = ሾߙ<sup>ଵ</sup>

0.3

0.4

0.5

0.6

0.7

0.8

Unit discharge, 11

0.9

1.0

1.1

thickness, ellipse ratio at the leading edge, and trailing edge. The design parameters for the stay vane are shown in Figure 3a. The design parameters for the stay vane are defined as in Equation (4).

50 55 60 65 70 75 80 85 90 95 100 105

GVO 11

GVO 14

GVO 38

GVO 41

GVO 17

GVO 20

GVO 23 GVO 26 GVO 29 GVO 32 GVO 35

GVO 8

Unit speed, <sup>11</sup>

92 93

91

90 89

<sup>87</sup> <sup>86</sup>

81

85

83 82

84

88

$$\mathbf{d}\_{SV} = \begin{bmatrix} \alpha\_{1\prime} & \dots \end{bmatrix} \begin{bmatrix} \delta\_{1\prime} \dots \delta\_{5\prime} \delta\_{1\prime} \dots \delta\_{5\prime} a\_{le\prime} a\_{le} \end{bmatrix}^T \tag{4}$$

where *dSV* is the design variables matrix of the stay vane, α*<sup>i</sup>* is the vane angle at *i*th section of the stay vane, δ*<sup>i</sup>* is the thickness at *i*th section of the stay vane, *ale* is the ellipse ratio at the leading edge (LE), *ate* is the ellipse ratio at the trailing edge (TE) and superscript *T* indicates the transpose of the design variables matrix. ௌ *α δ*

**Figure 3.** (**a**) Parametric design schematic view and (**b**) initial stay vane shape.

Figure 3b shows the initial stay vane shape using the vane angle and thickness distributions. The inlet vane angle is 30◦ , and the outlet vane angle is 33◦ . The thickness of the stay vane is 46 mm at the LE and 28 mm at the TE. The maximum thickness of the stay vane is 150 mm at 0.4 normalized distance from the LE. The vane angle and thickness distribution are the same throughout the stay vane.

#### 2.1.4. Casing Design

The spiral casing shape is dependent on the cross-section radii [20]. Figure 4a,b indicates the parametric design parameters and initial shape design of the casing, respectively. The parametric design of the casing shape is defined as in Equation (5).

$$\mathbf{d\_{CA}} = \begin{bmatrix} r\_0, r\_2, \dots, r\_{11} \end{bmatrix}^T \tag{5}$$

 where *dCA* is the design variables matrix of the casing, *r<sup>i</sup>* is the cross-section radius at *i*th section of the spiral casing.

Vane angle, *α*

 (°)

ݎሾ =

<sup>ଶ</sup>ݎ , ଵଵሿݎ,⋯, ்

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Normalized distance from to

Vane Angle Thickness

Thickness, *δ* (mm)

**Figure 4.** (**a**) Parametric design and (**b**) initial cross-section radius of casing (red line indicates measuring location). 0.0 0 30 60 90 120 150 180 210 240 270 300 330 360 Central angle, *θ* (°)

#### 2.1.5. Draft Tube Design with J-Groove Installation

*θ* The diffuser angle of the draft tube is generally determined in the range of 3◦ to 10◦ for minimum energy loss, which is an important design parameter for the discharge pressure recovery and flow uniformity in the draft tube [25]. The conceptual design of the draft tube shows that the diffuser angle of 3.5◦ provides maximum pressure recovery. Figure 5 shows the technical design of the draft tube. Moreover, in the case of the off-design condition, there exists swirl flow in the draft tube of the Francis hydro turbine, and J-Grooves can be an effective countermeasure of the flow instability in the draft tube [10]. J-Grooves are the grooves that are installed on the draft tube inner wall of the Francis turbine. The design parameters of the J-Groove are defined as angle (θ*JG*), length (*lJG*), depth (*dJG*) and number (*nJG*), which are shown in Figure 6. The J-Groove is used to suppress the swirl flow in the draft tube by the reverse flow mechanism through the J-Groove passage [11]. *θ*

**Figure 5.** *Cont*.

**Figure 5.** (**a**) Top view and (**b**) side view of draft tube (all dimensions are in mm).

൧ீ݊ ீ,݈ ீ,ߠ ீ,ൣ݀ = ் ் , **Figure 6.** J-Groove shape install on the draft tube inner wall and design parameters.

*θ* ் The parametric design for the J-Groove is represented in Equation (6).

$$\mathbf{d}\_{\rm DT} = \left[ d\_{\rm IG}, \theta\_{\rm IG}, l\_{\rm IG}, n\_{\rm IG} \right]^T,\tag{6}$$

ீ݊ ߠ ݊ீ = 180° *θ* ் ݊ீ where *dDT* is the design variable matrix of the J-Groove, *dJG* is the J-Groove depth, θ*JG* is the J-Groove angle, *lJG* is the J-Groove length and *nJG* is the number of J-Grooves.

ீߠ ீ݊ ߠThe grooves are evenly distributed in the circumference of the draft tube. Therefore, the number of J-Grooves can be calculated as in Equation (7). The initial shape of the J-Groove is defined as *dJG* = 106 mm, *lJG* = 2000 mm, θ*JG* = 12◦ and *nJG*= 15.

$$m\_{\rm fG} = \frac{180^{\circ}}{\theta\_{\rm fG}} \tag{7}$$

#### *2.2. Optimization Methodology*

#### 2.2.1. Process of Shape Optimization

<sup>ଵ</sup>ݔሺ݂ = ݕ <sup>ଶ</sup>ݔ , ݔ,⋯, ሻ ݁ The optimization process for the shape optimization of Francis hydro turbine fixed flow passages is illustrated in Figure 7. The design of experiments (DOE) was generated by using the optimal-space filling (OSF) method. The optimization for the fixed flow passages was carried out by using response surface methodology (RSM) and multi-objective genetic algorithm (MOGA). RSM is considered a sensitivity analysis tool, which is used to improve the sensitivity between the objective functions and input parameters [26]. RSM can be expressed as in Equation (8). The RSM uses a first-order and second-order polynomial form to develop the precise and concise correlation model, and the mathematical expressions of RSM are shown in Equations (9) and (10):

$$y = f(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n) + e \tag{8}$$

$$y = \alpha + \sum\_{i=1}^{n} a\_i \mathbf{x}\_i + e\_1 \tag{9}$$

$$y = a + \sum\_{i=1}^{n} a\_i \mathbf{x}\_i + \sum\_{i=1}^{n} a\_{ii} \mathbf{x}\_i^2 + \sum\_{i=1}^{n-1} \sum\_{j=i+1}^{n} a\_{ij} \mathbf{x}\_i \mathbf{x}\_j + e\_2 \tag{10}$$

݁<sup>ଵ</sup>

where *y* is the response of the system, *x*1, *x*2, · · · , *x<sup>n</sup>* are the independent variables and *e* is the error, *xi* is the *i*th input parameter, α is the coefficient of the response surface. ⋯ *α*

ݔߙ <sup>ఖ</sup>ߙ=ݕ 

ୀଵ

**Figure 7.** Optimization workflow for the 100 MW class Francis hydro turbine fixed flow passages. (DOE is Design of experiment and MOGA is Multi-Objective Genetic Algorithm).

RSM was used to decrease the computational cost in the optimization process. The response surface for the objective functions can be generated by using various RSMs such as genetic aggregation (GA) [27], radial basis function (RBF) [28], polynomial response surface (PRS) [29], Kriging (KG) [30], non-parametric regression (NPR) [31], neural network (NN) [32]. All these methodologies have their pros and cons. The selection of the RSM is based on the accuracy and consistency of the methodology. Among different RSMs, the accuracy is measured by using the goodness of fit test. The goodness of fit test can be calculated by using the coefficient of determination (CoD), maximum relative residual (MRR) and root mean square error (RMSE), which are expressed in Equations (11)–(13).

$$\text{CoD} = 1 - \frac{\sum\_{i=1}^{n\_s} (y\_i - \hat{y}\_i)^2}{\sum\_{i=1}^{n\_s} (y\_i - \overline{y}\_i)^2} \tag{11}$$

ሻ <sup>ೞ</sup> ଶ

$$MRR = \max\_{i} \left[ Abs\left(\frac{y\_i - \hat{y}\_i}{\overline{y}}\right)\right] \tag{12}$$

$$RMSE = \sqrt{\frac{1}{n\_s} \sum\_{i=1}^{n\_s} (y\_i - \mathfrak{g}\_i)^2} \tag{13}$$

where *n*<sup>s</sup> is the number of verification points, *y<sup>i</sup>* is the response from CFD analyses, *y*ˆ*<sup>i</sup>* is the corresponding response from the surrogate model and *y* is the arithmetic mean of *y*<sup>i</sup> . The verification points are used to evaluate Equations (11)–(13). If the result of the goodness measure shows CoD = 100%, MRR = 0% and RSME = 0%, it means that the response surface is highly precise and accurate.

The goodness measures concluded that the genetic aggregation was most suitable for the approximation of response surface, compared to other RSMs. The results of the goodness of fit test are shown in Table 2. Therefore, in this study, the genetic aggregation method was implied for the preparation of the response surface. The optimization was carried out using MOGA. Table 3 indicates the setting criteria for the optimization process.


**Table 2.** Results of goodness of fit test. (CoD is Coefficient of Determination, MRR is Maximum Relative Residual, RMSE is Root Mean Square Error).

**Table 3.** Information of setting criteria for MOGA.


#### 2.2.2. Process of Stay Vane Shape Optimization

In this study, the turbine efficiency η(*dSV*), flow uniformity γ(*dSV*), head loss *H<sup>l</sup>* (*dSV*), effective head *H*(*dSV*), and flow rate *Q*(*dSV*) were considered for the evaluation of stay vane design as in Equations (14)–(18). The measurement locations of the flow uniformity at *SVout* and head loss are calculated by the difference between the total pressure at *SVin* and *SVout*, as shown in Figure 3.

$$
\eta(\mathbf{d}\_{SV}) = \begin{bmatrix} \tau \omega \\ \rho \mathbf{g} \mathbf{Q} \mathbf{H} \end{bmatrix} \times 100\% \tag{14}
$$

$$\gamma(d\_{SV}) = \left[1 - \oint \frac{\sqrt{\left(\overline{u} - u\right)^2}}{2A\overline{u}} dA\right] \times 100\% \tag{15}$$

$$H\_l(\mathbf{d}\_{SV}) = \frac{\Delta p\_{total \oplus SV}}{\rho \mathbf{g}} \tag{16}$$

$$H(\mathbf{d}\_{SV}) = \frac{p^t\_{\text{inlet}} - p^t\_{\text{outlet}}}{\rho \mathbf{g}} \tag{17}$$

$$Q(d\_{SV}) = \frac{m\_{outlet}}{\rho} \tag{18}$$

where τ is torque generated by runner (Nm) and ω is rotational speed of runner (rad/s). *u* is the average velocity in stay vane passage (m/s), *u* is the local velocity in stay vane passage (m/s) and *A* is the cross-section area of stay vane passage (m<sup>2</sup> ). ∆*ptotal*@*SV* is the change in total pressure in stay vane passage (Pa), *p t inlet* and *p t outlet* are total pressures at inlet and outlet of the turbine (Pa). *moutlet* is mass flow rate of water at the turbine outlet (kg/s).

The optimization formulation for the stay vane is elaborated as in Equation (19).

$$\begin{array}{c} \text{maximize } \eta(\boldsymbol{d}\_{SV}), \ \boldsymbol{\gamma}(\boldsymbol{d}\_{SV})\\ \text{minimize } H\_l(\boldsymbol{d}\_{SV})\\ \text{subject to } 80 \text{ m} \le H(\boldsymbol{d}\_{SV}) \le 95 \text{ m} \\ 120 \text{ m}^3/\text{s} \le Q(\boldsymbol{d}\_{SV}) \le 135 \text{ m}^3/\text{s} \\ \boldsymbol{d}\_{SV}^\text{L} \le \boldsymbol{d}\_{SV} \le \boldsymbol{d}\_{SV}^\text{U} \end{array} \tag{19}$$

where *d* L *SV* and *d* U *SV* are lower and upper bounds for the design variable *dSV*, respectively, and their values are summarized in Table 4. The turbine efficiency η(*dSV*) and flow uniformity γ(*dSV*) were maximized to obtain more output power, while vortex-induced vibration was suppressed. At the same time, the head loss *H<sup>l</sup>* (*dSV*) was minimized to prevent loss of power in the stay vane. The effective head *H*(*dSV*) and flow rate *Q*(*dSV*) of the turbine were used as constraints for the stay vane design.


**Table 4.** Bounds for design variables of stay vane.

The optimal Pareto front for the stay vane design is shown in Figure 8. The Pareto front is plotted between turbine efficiency and flow uniformity. The trade-off between turbine efficiency and flow uniformity is required to obtain the optimal design of stay vane. Flow uniformity is measured at the outlet of the stay vane. The measuring location plays a vital role in the calculation of flow uniformity. If the measuring location changes, the nature of the Pareto front will change. The main objective of the design optimization is to have smooth flow distribution in the stay vane flow passage. The flow uniformity of the optimal stay vane should be above 90%. Based on the above assumptions, the optimal stay vane (OSV) was selected with flow uniformity 91.97% and turbine efficiency 96.37%.

**Figure 8.** Pareto front for the optimization of stay vane design.

ଶ

Δ௧௧@௦ ݃ߩ

ሺሻ subject to ߛሺሻ ≥ 97%

തݑܣ2

maximize ߛሺሻ ܪ minimize

> ≥ ≥

100% × ൩ܣ݀

∆௧௧@௦

ߛሺ<sup>ሻ</sup> = 1 − ර ඥሺݑത−ݑሻ

ܪ ሺሻ =

ܪ ሺሻ

തݑ

  ሻሺߛ

#### 2.2.3. Process of Casing Shape Optimization

In order to evaluate the flow condition in the casing, the flow uniformity γ(*dCA*) and head loss *Hl* (*dCA*) were examined. The flow uniformity, which determines the deviation of flow velocity in the casing shape, is calculated as in Equation (20). The flow uniformity was measured as an averaged value at the location of the whole casing outlet of 1.00 *D<sup>e</sup>* from the runner axis center, which is shown in Figure 4a by a red circle. The head loss was defined by the losses in the spiral casing passage due to flow mixing and wall friction, as in Equation (21), and the head loss was calculated by the difference between inlet and outlet of casing.

$$\gamma(d\_{CA}) = \left[1 - \oint \frac{\sqrt{(\overline{u} - u)^2}}{2A\overline{u}} dA \right] \times 100\% \tag{20}$$

$$H\_l(\mathbf{d\_{CA}}) = \frac{\Delta p\_{total\\$cooling}}{\rho \mathbf{g}} \tag{21}$$

where *u* is the average velocity in casing passage (m/s), *u* is the local velocity in casing passage (m/s), and *A* is the cross-section area of casing passage (m<sup>2</sup> ), ∆*ptotal*@*casing* is change in total pressure in casing passage (Pa).

The design optimization problem of the casing is formulated as in Equation (22).

$$\begin{array}{l}\text{maximize } \gamma(\mathsf{d}\_{\mathsf{CA}})\\\text{minimize } H\_{l}(\mathsf{d}\_{\mathsf{CA}})\\\text{subject to } \gamma(\mathsf{d}\_{\mathsf{CA}}) \ge 97\%\\\ \mathsf{d}\_{\mathsf{CA}}^{\mathsf{L}} \le \mathsf{d}\_{\mathsf{CA}} \le \mathsf{d}\_{\mathsf{CA}}^{\mathsf{U}}\end{array} \tag{22}$$

where *d* L *CA* and *d* U *CA* are lower and upper bounds for the design variable *dCA*, respectively.

The bounds for design variables of the casing are shown in Table 5. The Pareto front for the casing shape optimization was prepared by using head loss and flow uniformity. Figure 9 shows the Pareto front for the optimization of the casing. The Pareto front shows the trade-off between flow uniformity and head loss. Thus, the selection of the optimal design is based on the requirement of the user. In this study, the main objective was to increase the flow uniformity above 97%. Therefore, the optimal design was selected considering flow uniformity above 97% with minimum head loss.


൫ ۺ

൯ ൫

܃ ൯

**Figure 9.** Pareto front for the optimization results of casing design.

#### 2.2.4. Process of Draft Tube Shape Optimization

ܪ ሻ்ሺܵ ሻ்ሺߟ ሺ்ሻ ܪሺ்ሻ ܳሺ்ሻ The optimization was carried out to obtain the optimal solution for the draft tube shape with the J-Groove installation. The turbine efficiency η(*dDT*), swirl intensity *S*(*dDT*), head loss *H<sup>l</sup>* (*dDT*), effective head *H*(*dDT*), and flow rate *Q*(*dDT*) were considered for the evaluation of draft tube shape design, as in Equations (23)–(25). The measurement locations of the swirl intensity were set in the range of *z*/*R*<sup>0</sup> = 1.15 to 3.60, as shown in Figure 6.

$$
\eta(\mathbf{d}\_{\rm DT}) = \left[\frac{\tau \omega}{\rho g \mathbf{Q} H}\right] \times 100\% \tag{23}
$$

$$S(d\_{DT}) = \frac{\int v\_{\theta} v\_{a} r^{2} dr}{R \int v\_{a}^{2} r dr} \tag{24}$$

$$H\_l(\mathbf{d}\_{\rm DT}) = \frac{\Delta p\_{\rm total\#IG}}{\rho g} \tag{25}$$

*τ* ω *ρ* ℃ ݒ <sup>ఏ</sup>ݒ ∆௧௧@ீ where τ is torque generated by runner (Nm), ω is rotational speed of runner (rad/s), ρ =997 kg/m<sup>3</sup> is the density of water at 25 °C, *g* = 9.81 m/s 2 is gravitational acceleration, *Q* is the flow rate (m<sup>3</sup> /s), and *H* is the effective head (m), *v*<sup>θ</sup> is the local tangential velocity in the draft tube (m/s), *v<sup>a</sup>* is the local axial velocity in the draft tube, *r* is the radial position, *R* is the cross-section radius, ∆*ptotal*@*JG* is the change in total pressure in draft tube passage (Pa), *z* is the vertical distance from the axis of turbine, *R*<sup>0</sup> is the runner outlet radius.

ܵሺ்ሻ In order to investigate the flow instability and to express the complicated and unique internal flow behavior in the draft tube effectively, swirl intensity *S*(*dDT*) was adopted to determine the strength of swirl flow in the draft tube. The swirl intensity represents the ratio of the axial flux of angular momentum to axial momentum, as shown in Equation (24).

The optimization formulation for the draft tube is elaborated as in Equation (26).

$$\begin{array}{c} \text{maximize } \eta(\mathbf{d}\_{\text{DT}})\\ \text{minimize } \mathcal{S}(\mathbf{d}\_{\text{DT}}), \ H\_{l}(\mathbf{d}\_{\text{DT}})\\ \text{subject to } 80 \text{ m} \le H(\mathbf{d}\_{\text{DT}}) \le 95 \text{ m} \\\ 120 \text{ m}^3/\text{s} \le \mathcal{Q}(\mathbf{d}\_{\text{DT}}) \le 135 \text{ m}^3/\text{s} \\\ \mathbf{d}\_{\text{DT}}^\text{L} \le \mathbf{d}\_{\text{DT}} \le \mathbf{d}\_{\text{DT}}^\text{U} \end{array} \tag{26}$$

where *d* L *DT* and *d* U *DT* are lower and upper bounds for the design variable *dDT*, respectively, and their values are summarized in Table 6. The turbine efficiency η(*dDT*) was maximized to obtain more output power. At the same time, the swirl intensity *S*(*dDT*) and head loss *H<sup>l</sup>* (*dDT*) were minimized to suppress ்

்

݊ீ

the swirl flow and prevent energy loss in the draft tube flow passage. The effective head *H*(*dDT*) and flow rate *Q*(*dDT*) of the turbine were used as constraints for the draft tube design, which are expressed as in Equations (17) and (18), respectively. ሻ்ሺߟ ܪ ሻ்ሺܵ ሺ்ሻ ܪሺ்ሻ ܳሺ்ሻ

்

maximize ߟሺ்ሻ minimize ܵሺ்ሻ, ܪ

subject to 80 m ≤ ܪሺ்ሻ ≤ 95 m 120 m<sup>ଷ</sup>⁄s ≤ ܳሺ்ሻ ≤ 135 m<sup>ଷ</sup>⁄s

் ≥ ் ≥

்

ሺ்ሻ

,


**Table 6.** Bounds for design variables of draft tube shape.

The lower and upper limits of draft tube design variables are indicated in Table 6.

The optimization of the draft tube design was carried out at the design point. The Pareto front was prepared by the trade-off between turbine efficiency and swirl intensity, which is shown in Figure 10.

**Figure 10.** Pareto front for optimization of draft tube design at the design point.

#### **3. CFD Methodology**

The CFD analysis for the turbomachinery requires a highly reliable computational system for the calculation of complex internal flow phenomena. Moreover, while conducting the optimization, numerous samples are needed, which demand extensive computational cost and time for CFD analysis. Figure 11 shows the numerical scheme of CFD analysis adopted in this study, in combination with the optimal design process. The CFD analysis process is directly connected to the optimum design process. The CFD analysis method was adopted from previous studies [33–35].

The CFD analysis for the casing DOE samples was performed without stay vanes because the flow field in the spiral casing is independent of the flow field of the stay vane [36]. Furthermore, the single flow passage analysis for the stay vane was used for the calculation of DOE. It provides precise CFD analysis results and reduces the computation time. However, the full domain analysis was required for the DOE of the draft tube with the J-Groove installation because the flow field in the J-Groove is dependent on other components of the Francis hydro turbine.

The CFD analysis was conducted using a commercial code of ANSYS CFX 19.2 [37]. The numerical analysis was performed by solving the governing equations and Reynolds-averaged Navier–Stokes (RANS) with the turbulence model. In this study, the Shear Stress Transport (SST) turbulence model was selected because the SST model combines the capabilities of the κ*-*ω model away from the walls and the robustness of the κ*-*ε turbulence near the walls by using blending functions of the automatic near-wall treatment. The Rhie–Chow algorithm was used to interpolate the pressure–velocity coupling mechanism. The high-resolution order was used to solve the advection term, and the first-order upwind difference was used to solve the turbulence numeric [37].

*κ ε*

*κ ω*

**Figure 11.** Numerical scheme of CFD analysis in combination with optimal design process (RANS is Reynolds-averaged Navier-Stokes).

The proper numerical grids are required for precise and accurate computational analysis. The structured mesh for the numerical analysis was generated using ANSYS ICEM 19.2 [37]. The numerical grids for the 100 MW class Francis hydro turbine are shown in Figure 12. The mesh dependence test was carried out to determine the optimum number of nodes. The results of the mesh dependence were compared with efficiency and output power. Figure 13 shows the mesh dependence test results, and we concluded that 8.2 million nodes was the optimum number for computational analysis. Table 7 presents the information on the numerical grids used for CFD analysis. The non-dimensional wall distance *y*+ values for the several components of the Francis hydro turbine were less than 100, which was suitable for the SST turbulence model with automatic near-wall treatment within the reliable resolution range of 1 < *y* + < 100. Table 8 shows the summary of boundary conditions for CFD analysis. The performance curves of the 100 MW Francis hydro turbine are shown in Figure 14. The performance curves indicated that the design point and best efficiency point were matched well. They verified that the conceptual design of the 100 MW Francis hydro turbine was acceptable.

**Figure 12.** Numerical grids of the 100 MW class Francis turbine.

5 6 7 8 9 10 11 12 13 14

Number of nodes (in millions)

60

70

80

Efficiency, *η* (%)

90

100

−

80

85

Efficiency Power

90

Power, (MW)

95

100

Casing

Guide Vane

Stay Vane

Runner

Draft Tube

Thickness, *δ*

(mm)

**Figure 13.** Mesh dependency test results for CFD analysis.


**Table 8.** Summary of boundary conditions for CFD analysis.


**Figure 14.** Validation of conceptual design of the 100 MW Francis hydro turbine by CFD analysis.

Vane Angle [Initial Stay Vane] Vane Angle [Optimal Stay Vane] Thickness [Initial Vane Angle] Thickness [Optimal Stay Vane]

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Initial Stay Vane [ ] Optimal Stay Vane [ ]

Normalized distance from LE to TE

Vane angle,

*α* (°)

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

#### *4.1. Stay Vane Shape Optimization*

0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%

Efficiency, *η*

Figure 15 shows the comparison of the vane angle and thickness between the initial stay vane (ISV) and optimal stay vane (OSV) shapes. The inlet vane angle was changed from 30◦ to 32◦ . The position for the maximum thickness was modified from the normalized distance of 0.4 to 0.3. The maximum thickness was increased from 150 mm to 158 mm. The 3D view of the initial and optimal stay vane shapes are shown in Figure 15. Table 9 shows the results of stay vane optimization for the 100 MW class Francis turbine. The targeted objectives of turbine efficiency, flow uniformity and head loss were all improved remarkably.

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2

Design Point

Efficiency Power

Unit Discharge, <sup>11</sup>

Power, (MW)

**Figure 15.** Comparison of initial and optimal stay vane designs for the 100 MW class Francis turbine.


**Table 9.** Results of stay vane optimization for the 100 MW class Francis turbine.

The flow uniformity was measured at the outlet of the stay vane. The flow uniformity encountered the average deviation of local velocity in the reference area. The flow deviation in the small section did not show a significant effect on flow uniformity. Therefore, flow angle distribution and vorticity were evaluated for internal flow patterns in stay vane passage. The comparison of the flow angle between the ISV and OSV is shown in Figure 16. The flow angle (θ*u*) is defined by Equation (27).

$$\theta\_{\rm \ell} = \tan^{-1} \left( \frac{v\_{\ell}}{v\_r} \right) \tag{27}$$

where *v*<sup>θ</sup> and *v<sup>r</sup>* are tangential and radial velocity components.

Figure 16 shows the comparison of flow angle in the ISV and OSV at the design point. The peak value of the flow angle at the hub and shroud indicated the occurrence of the secondary flow near the hub and shroud walls. The flow angle difference in the passage between the outlet of stay vane and inlet of guide vane in the OSV became remarkably smaller than that of the ISV. It meant that the OSV had a relatively larger ability to maintain a proper flow angle. The smaller difference caused the lower vorticity in the passage between the stay vane and guide vane in the OSV.

*θ*

ߠ<sup>௨</sup> = tanିଵ ൬

ఏݒ ݒ ൰ *θ*

**Figure 16.** Comparison of flow angle in the ISV and OSV at the design point.

Figure 17 shows the strength of vorticity in between the cascade passages of the stay vane and guide vane. The decrease in the vorticity at the OSV flow passage made the flow smoother. Thus, the possibility of occurrence of secondary flow and vortices at the OSV flow passage decreased significantly in the OSV. Therefore, the OSV made a more uniform flow distribution in the vane's passage.

**Figure 17.** Comparison of vorticity in ISV and OSV flow passages (**a**) 0.25 span and (**b**) 0.75 span at design point.

#### *4.2. Casing Shape Optimization*

0

500

1000

1500

Cross-section radius, (mm)

2000

2500

*θ θ* Figure 18 indicates the cross-section radius comparison between the initial and optimal shape of the casing. The cross-section radius of the optimal casing shape was greater than the initial casing shape at the central angles below θ = 180◦ , but the cross-section radius near the casing tongue of θ = 345◦ was almost the same.

Initial Design Optimal Design

0 30 60 90 120 150 180 210 240 270 300 330 360

Central angle, *θ* (°)

*θ θ*

**Figure 18.** Comparison of initial and optimal casing design for the 100 MW class Francis turbine.

Table 10 shows the flow uniformity and head loss by the initial and optimal casing shapes, which was compared to the design point. The flow uniformity increased slightly in the optimal casing in comparison with that of the initial casing shape; furthermore, the optimal casing design showed a significant decrease in the head loss in comparison with that of the initial casing design.

**Table 10.** Results of casing optimization for the 100 MW class Francis turbine.


The secondary vortex intensity *J n ABS* was used to evaluate the internal flow behavior in the casing quantitatively. The area integral of vorticity around the casing was calculated as defined in Equation (28).

$$J\_{ABS}^n = \frac{1}{A} \iint\limits\_0 \left| \nabla \times \overrightarrow{u}^n \right| dA \tag{28}$$

where <sup>→</sup> *u n* is the velocity at normal direction to the cross section, *A* is the cross section area.

Figure 19 shows the comparison of secondary vortex intensity between the initial and optimal casing designs. The secondary vortex intensity is in increasing order from the inlet to the casing tongue. The vortex intensity was suppressed significantly by optimal design in comparison to the initial casing shape. Therefore, it was clear that the flow uniformity and head loss could be improved effectively by the current optimum design process.

#### *4.3. Draft Tube Shape Optimization*

Table 11 shows the comparison of the design parameter sizes of the J-Grooves for draft tube shape optimization. The 3D model of the J-Grooves installed on the draft tube wall is shown in Figure 20.

**Table 11.** Comparison of design parameter size of J-Grooves for draft tube shape optimization.


0.50

0.55

Ԧሬݑ 

ௌܬ 

> ௌܬ = 1 ܣ

ඵ |∇×ݑሬԦ


ܣ݀

**Figure 19.** Comparison of secondary vortex intensity between initial and optimal casing designs.

*θ* **Figure 20.** Comparison between (**a**) initial (**b**) optimal J-Groove models installed on the draft tube walls.

݊ீ Figure 21 shows the comparison of the swirl intensity in the draft tube by J-Groove shapes at the design point. The swirl intensity was suppressed significantly with the installation of J-Groove. There was a 12.12% swirl intensity reduction with initial J-Groove installation from the case without J-Groove installation. Moreover, the additional 6.64% swirl intensity reduction was achieved by the optimized J-Groove shape from the initial J-Groove shape. Therefore, it was clear that the installation of an optimal J-Groove in the draft tube had a significant effect on the suppression of the flow instability caused by the swirl flow. Table 12 reveals the results summary of the draft tube shape optimization. The study results indicated that the installation of the J-Groove on the wall of the draft tube had almost no influence on the turbine performance but suppressed the flow instability of the swirl flow remarkably.

**Table 12.** Summary of draft tube shape optimization results for the 100 MW class Francis turbine.


1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8

0

Without J-Groove With Initial J-Groove With Optimal J-Groove

**Figure 21.** Comparison of swirl intensity in the draft tube by J-Groove shapes at the design point.

#### **5. Conclusions**

In the present study, the fixed flow passage shapes of a 100 MW class Francis hydro turbine were optimized for the internal flow uniformity by a CFD-based shape design optimization process. The stay vane, casing, and draft tube were optimized separately to understand the flow characteristics in each flow passage. The objective of the optimization was to maximize the flow uniformity and minimize the head loss in each flow passage.

A CFD-based shape design optimization process of the parametric conceptual design, detailed design, and optimal design of the fixed flow passage of the Francis hydro turbine was accomplished. The design and optimization process can be generalized for the reaction hydro turbine stay vane, casing, and draft tube with J-Grooves. Moreover, better flow uniformity was achieved in the Francis hydro turbine by the fixed flow passages optimization process. For the optimization process, response surface methodology was used to generate the response surface, and a multi-objective genetic aggregation method was used to determine the global optimum solution via the optimal Pareto front.

The optimum stay vane shape was achieved with the remarkably decreased vorticity around the stay vane flow passage, which resulted in the highly improved flow uniformity in the vane passage. The optimal casing passage shape was achieved with the increased flow uniformity and the significantly decreased head loss in comparison with that of the initial casing shape. The secondary vortex intensity was suppressed effectively by casing shape optimization. The installation of a J-Groove on the wall of the draft tube had almost no influence on the turbine performance but suppressed the flow instability of swirl flow remarkably in the draft tube passage.

**Author Contributions:** Conceptualization, U.S., and Y.-D.C.; methodology, U.S.; software, Y.-D.C.; validation, U.S. and Y.-D.C.; formal analysis, U.S.; investigation, U.S.; resources, Y.-D.C.; data curation, U.S.; writing—original draft preparation, U.S.; writing—review and editing, U.S. and Y.-D.C.; visualization, U.S.; supervision, Y.-D.C.; project administration, Y.-D.C.; funding acquisition, Y.-D.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the New and Renewable Energy of the Korea Institute of Energy Technology Evaluation and Planning (KETEP), which was funded by the Korean Government Ministry of Trade, Industry and Energy (no. 20163010060350).

**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* **Aerodynamic Performance of an Octorotor SUAV with Di**ff**erent Rotor Spacing in Hover**

### **Yao Lei 1,2,\* , Yuhui Huang <sup>1</sup> and Hengda Wang <sup>1</sup>**


Received: 3 October 2020; Accepted: 26 October 2020; Published: 28 October 2020

**Abstract:** To study the aerodynamic performance of hovering octorotor small unmanned aerial vehicles (SUAV) with different rotor spacing, the computational fluid dynamics (CFD) method is applied to analyze the flow field of an octorotor SUAV in detail. In addition, an experimental platform is built to measure the thrust and power of the rotors with rotor spacing ratios *L*/*D* of 1.0, 1.2, 1.4, 1.6, and 1.8, sequentially. According to the theory of momentum, rotor aerodynamic performance is obtained with qualitative analysis. Further analysis with numerical simulation is presented with the flow field of the octorotor SUAV, the vorticity distribution, velocity distribution, pressure distribution, and streamline. The results show that the aerodynamic performance varies with the rotor spacing. Specifically, the aerodynamic performance is poor at *L*/*D* = 1.0, which is accompanied with strong interaction of wake and tip vortexes and interaction with each other. However, the aerodynamic efficiency is much improved with a larger rotor spacing, especially achieving the highest at *L*/*D* = 1.8, which is considered to be the best rotor spacing ratio for this kind of octorotor SUAV.

**Keywords:** octorotor SUAV; aerodynamic performance; rotor spacing; hover; CFD; vortices distribution

### **1. Introduction**

Small unmanned aerial vehicles (SUAVs) are normally less than 25 kg and easy to pack. Especially, the octorotor SUAVs with evenly distributed rotors have been widely used in agricultural, surveillance, and military, due to its advantages of simple operation and convenient portability. In addition, octorotor SUAVs have higher load and more damaged redundancy as compared with a quadrotor or hex-rotor SUAVs.

Considering that octorotor SUAVs often operate in environments where the Reynolds numbers are less than 10<sup>5</sup> , viscosity, laminar separation bubbles, thickened boundary layer, and flow separation influence the rotor tip, which can lead to increased drag on the vehicle, thus, resulting in poor aerodynamic performance. Especially, there is a strong vibration at higher rotor speed which can affect the dynamic stability of a rotorcraft [1–3]. There are two main ways to improve the aerodynamic performance of small rotary-wing UAVs. One way is to change the parameters of the blade [4], such as changing the shape of the blade platform, adding a twist, taper, and so on. The second way is to change the layout and configuration of the rotors, such as changing the distance between adjacent rotors, the number of rotors, the tilt angle of the rotors, and so on. At the same time, with the help of constantly developing computer simulation technology, we can intuitively explore the reasons that lead to a decrease or increase in aerodynamic performance of the rotor, and help to optimize the layout of the SUAVs.

Early studies have focused on hover performance analysis of a small single rotor and a small coaxial rotor. Bohorquez [5–7] experimentally measured the aerodynamic performance of the rotor and used fluorescent oil to visualize surface flow, which showed consistently poor performance for a variety of rotors with low Reynolds numbers (less than Re*tip* <sup>=</sup> 0.5 <sup>×</sup> <sup>10</sup><sup>5</sup> ). Ramasamy [8] used the digital particle image velocimetry (DPIV) to study the aerodynamic efficiency of the hovering micro rotor with different blade shapes and found that changing the shape of the blade, such as adding a linear twist and changing the planform of the airfoil, could effectively improve the aerodynamic performance of the blade. A series of parametric studies of hovering coaxial rotor were conducted by Lakshminarayanan [9] and Syal [10]. Changing the inter-rotor spacing proved to influence the aerodynamic of the coaxial rotor in hover.

In recent years, research on aerodynamic performance of small rotorcrafts has deepened. Yoons [11] studied the influence of the turbulence model on the flow solution accuracy of hovering rotor. High-fidelity computational fluid dynamics (CFD) simulations have been implemented for multi rotor [12,13], and adding components under the airframe was found to weaken the interactions between rotors. Henricks [14,15] investigated the effect of rotor variables on the aeroacoustics performance of a hovering rotor. Greater twist and taper were proven to improve both the aerodynamic and acoustic performance. A high-speed stereo particle-image velocimetry (SPIV) study on multi rotor was conducted by Shukla [16,17], and vortex–vortex, blade–vortex, and vortex–duct interactions were visualized. The effect of a tilting rotor on the aerodynamic performance was studied by Zhang [18]. Tilting has been shown to help keep the vehicle stable. The simulation of the downwash flow of octorotor was done by Yang [19].

The current study of rotor aerodynamic performance concentrated on the coaxial rotor and quadrotor since the rotor interference of an octorotor SUAV is relatively complicated, especially the different vortex–vortex and blade–vortex interactions with different rotor spacing. Therefore, analysis of the aerodynamic performance of octorotor is challenging work, which can also provide significant guidance for designing octorotor SUAVs.

The main research objective of this paper is to explore how the aerodynamic performance of a small octorotor SUAV changes with a change of rotor spacing, why it changes, and whether there is an optimal rotor spacing. In comparison with previous studies on small quadrotor and small hex-rotor SUAVs [20–22], the improvement of this study lies in the addition of preliminary studies on wake vortices.

#### **2. Theoretical Analysis**

Momentum theory is assumed to be an effective method for flow-field analysis of rotorcrafts. It regards the rotating rotor as an actuator disk and relates the induced airflow velocity through the actuator disk to the thrust and power of the rotor. According to the momentum theory [23], the induced flow velocity can be expressed as:

$$V\_h = \sqrt{\Gamma / 2\rho A} \tag{1}$$

where *<sup>T</sup>* is the thrust generated by the rotor, <sup>ρ</sup> is the air density, with a value of 1.225 kg · <sup>m</sup>−<sup>3</sup> , and *A* is the area of the rotor disk.

The induced power consumption can be expressed as:

$$P\_l = TV\_h = T^{3/2} / \sqrt{2\rho A} \tag{2}$$

To measure the aerodynamic efficiency of a rotor in hover, *FM* is defined as the ratio of induced power consumption to actual power consumption:

$$\text{F}M = \frac{P\_i}{P} = \frac{T^{3/2} / \sqrt{2\rho A}}{P} = \frac{C\_T^{3/2} / \sqrt{2}}{C\_P} \tag{3}$$

where

$$\mathcal{C}\_T = \frac{T}{\rho A \Omega^2 R^2} \tag{4}$$

is the thrust coefficient,

$$\mathcal{C}\_P = \frac{P}{\rho A \Omega^3 R^3} \tag{5}$$

is the power coefficient. Ω is the angular velocity of the rotor and *R* is the radius of the rotor.

*FM* of a given rotor is often expressed as a function of *CT*/σ. Here, σ is the solidity of the rotor, and can be defined as: 3 3

$$
\sigma = \text{Nc} / \pi \text{R} \tag{6}
$$

/

where *N* is the number of blades and *c* is the chord.

Since the tapered blade is adopted in this paper, the influence of varying chord length on the solidity must be considered, and Equation (7) is adopted to calculate the equivalent solidity: /

2

$$
\sigma\_{\ell} = 3 \int\_0^1 \sigma r^2 dr \tag{7}
$$

The total efficiency of the aircraft is determined by the power loading: 1

$$PL = \frac{T}{P} = \frac{\mathcal{C}\_T}{\mathcal{C}\_P \Omega R} \tag{8}$$

Combined with momentum theory, *PL* can also be expressed in the following form: =

$$PL = \frac{\sqrt{2\rho}FM}{\sqrt{DL}}\tag{9}$$

where *DL* is the disk loading and can be defined as: =

$$\text{DL} = \text{T} / \text{A} \tag{10}$$

In addition, the tip chord Reynolds number in this paper is defined by the chord length at 0.75*R*: = /

$$\mathbf{Re}\_{tip} = \rho V \mathbf{c}\_{0.75} / \mu \tag{11}$$

where <sup>µ</sup> is kinematic viscosity coefficient of air, with a value of 1.79 <sup>×</sup> <sup>10</sup>−<sup>5</sup> Pa · <sup>s</sup> with a temperature of 293.15 K, and a pressure of 101.325 kPa. 5 1.79 10

Figure 1 shows the arrangement of octorotor and interference between adjacent rotors.

**Figure 1.** Sketch of the octorotor arrangement.

#### **3. Hover Experiment**

#### *3.1. Experimental Setup*

The rotor used in this paper is made of carbon fiber, and unidirectional carbon fiber cloth is used as reinforcement. The weight of one rotor is 15 g. The rotor is operating at Re*tip* <sup>=</sup> 0.59 <sup>×</sup> <sup>10</sup>5~0.99 <sup>×</sup> 10<sup>5</sup> . The specific rotor parameters are shown in Table 1.



The schematic diagram of the test bench is shown in Figure 2. Each rotor is fitted with a separate motor and thrust sensors (model CZL601, accuracy ±0.02% F.S., range 0~3 kg). Adjacent rotors rotate in opposite directions, and the rotor speed is measured with an optical tachometer (model DT2234C, accuracy ±0.05 n% + 1 d). To reduce the impact brought by the ground effect, the rotor is installed at a height of 1.5 m above the ground. Rotor spacing *L* is defined as the distance between the rotor centers of two adjacent rotors. In this experiment, *L* = 400, 480, 560, 640, and 720 mm, respectively. Thus, the corresponding dimensionless spacing ratios *L*/*D* = 1.0, 1.2, 1.4, 1.6, and 1.8 (*D* is the diameter of the rotor). By changing the tip chord Reynolds number, the thrust and power consumption generated by the octorotor at five spacing are collected, and the aerodynamic performance of the octorotor is analyzed accordingly.

**Figure 2.** Schematic diagram of test bench.

#### *3.2. Error Analysis*

The main sources of error in the experiments are the standard deviations of the rotational speed and the mean voltages from the force sensors. Typical values of the standard deviations of thrust are about 1% of the mean values. Rotational speed error is related to the finite number of magnets that excite the tachometer, which causes an error of 1/24 × 60 = 2.5 RPM. The statistical (random) error in the coefficient values for each run was estimated by the standard deviation of repeated samples. This experiment mainly measured variable as thrust *T*, torque *Q*, and rotational speed Ω. To determine the accuracy of the test measurements, a normalized criterion is needed in this case [24–26]. Consider a

1

1/2 <sup>2</sup>

variable *X<sup>i</sup>* , the uncertainty of the calculated results can be expressed by using the "Kline–McClintock"' equation [27]:

$$
\Delta \mathcal{R} = \left\{ \sum\_{i=1}^{N} \left( \frac{\partial \mathcal{R}}{\partial X\_i} \Delta X\_i \right)^2 \right\}^{1/2} \tag{12}
$$

where ∆*X<sup>i</sup>* is the uncertainty of *X<sup>i</sup>* .

The expression of uncertainty of power can be obtained as follows:

$$
\Delta P = \sqrt{\left(\frac{\partial P}{\partial Q} \Delta Q\right)^2 + \left(\frac{\partial P}{\partial \Omega} \Delta \Omega\right)^2} \tag{13}
$$

The uncertainty of *P* as a percentage is:

$$\frac{\Delta P}{P} = \sqrt{\left(\frac{\Delta Q}{Q}\right)^2 + \left(\frac{\Delta \Omega}{\Omega}\right)^2} \tag{14}$$

Similarly:

$$\frac{\Delta PL}{PL} = \sqrt{\left(\frac{\Delta T}{T}\right)^2 + \left(-\frac{\Delta Q}{Q}\right)^2 + \left(-\frac{\Delta \Omega}{\Omega}\right)^2} \tag{15}$$

$$\frac{\Delta FM}{FM} = \sqrt{\left(\frac{3}{2}\frac{\Delta T}{T}\right)^2 + \left(-\frac{\Delta Q}{Q}\right)^2 + \left(-\frac{\Delta \Omega}{\Omega}\right)^2} \tag{16}$$

By calculation, the uncertainty of *P*, *PL,* and *FM* is 1.41%, 1.44%, and 1.51%, respectively. The values of uncertainty that are presented in this study are all calculated for 95% confidence levels.

#### *3.3. Results*

The thrust and power consumption of the octorotor with different spacing ratios is shown in Figure 3.

**Figure 3.** Thrust and power consumption variation.

It can be seen from Figure 3 that at the same power consumption, the thrust is significantly greater than the simple superposition of the thrust generated by the eight small single rotors. In addition, the thrust generated by the octorotor at *L*/*D* = 1.0 is significantly lower than that generated by the octorotor at the other four spacing ratios. It indicates that the rotor interference with proper rotor spacing is beneficial to improve the hover efficiency of the octorotor.

The variation of *FM* with different rotor spacing ratios is shown in Figure 4.

**Figure 4.** Variation of the figure of merit.

It can be found from Figure 4 that the *FM* of the octorotor shares a similar tendency, and it is very interesting to note that there is a sudden drop at *Re* = 79246 for *L*/*D* = 1.2 and at *Re* = 71,321 for *L*/*D* = 1.0. The likeliest explanation for this phenomenon is the collapse of the suction forces on the tip and the increase in vibration caused by the flow separation, which are relatively unsteady for smaller rotor spacing, where the rotor is apt to have somewhat greater interaction with its own wake. This heightened interaction is reflected in the greater variability to a decrease in *FM*. Furthermore, because of the unpredictable inflow and outflow of smaller rotor spacing where the rotor is apt to have interaction with its own wake, it is, therefore, not entirely surprising that some data do not follow the same trends. Additionally, the aerodynamic performance at *L*/*D* = 1.8 is much better than other spacing ratios with a much higher *FM*.

The variation of power loading versus disk loading is shown in Figure 5.

**Figure 5.** Variation of power loading and disk loading.

 (N/m<sup>2</sup> )

As can be seen from Figure 5, *PL* value at *L*/*D* = 1.8 always is higher than the other rotor spacings. The *PL* for 1.4 and 1.2 are approximately the same with a lower value. Especially, the *PL* for *L*/*D* = 1.0 was presented with lowest value which suffered from strong rotor interference.

To sum up, when the spacing ratio is 1.8, the small octorotor has better aerodynamic efficiency and total efficiency, and therefore it can be regarded as the best rotor spacing in the current working condition.

The *FM* comparison of the small octorotor, the small quadrotor [21], and the small hex-rotor [20] at the optimum rotor spacing is shown in Figure 6.

**Figure 6.** The *FM* comparison with quadrotor and hex-rotor.

It can be seen from Figure 6 that the aerodynamic performance of a small octorotor is much better than that of a small quadrotor and a small hex-rotor composed of the same rotor.

#### **4. Numerical Simulation**

#### *4.1. Simulation Setting and Mesh Analysis*

To further analyze the aerodynamic performance of the small octorotor, the flow field of the small octorotor with different spacing ratios was numerically simulated with ANSYS when Re*tip* <sup>=</sup> 0.79 <sup>×</sup> <sup>10</sup><sup>5</sup> . The detailed mesh distribution is shown in Figure 7.

**Figure 7.** Mesh distribution.

The whole computational domain is divided into nine regions including one cylinder stationary region and eight cylinder rotating regions (to capture the flow detail of four rotors with refined mesh), which has a total size of 70 million cells. In addition, the SUAV is located in the upper region of the domain to obtain the detail of the downwash flow of the SUAV. To handle the low Reynolds number flows, multiple viscous boundary layers are used in the simulations, and the max element metrics is below 0.8 to capture the flow detail of the rotor tip and the interfaces between stationary and rotating regions. The grid is considered to be sufficient to handle these simulations and the mesh on the rotor tip is refined to reach the independence state. For the boundary conditions, all three surfaces of the cylinder domain are set as no slip wall where the rotating region and stationary overlap is set as interface. The sliding mesh is used for the transient solver. The Spalart–Allmaras turbulence model was chosen for the RANS closure which was proven to be capable of handling the aerodynamic performance of the

multi-rotors, especially at low *Re*. The physical time step corresponds to a rotor rotation of 30 degrees, and 12 steps correspond to a rotor revolution. The rotor made a total of 40 revolutions. Furthermore, the mesh independence study is conducted to show that the results are already reaching the grid independence state.

For the numerical simulation, the accuracy of the model used, the quality of the mesh, the selected boundary conditions, the turbulence model, and the discrete format all produce some errors. For the experiment, the sensor error is the main error source. The comparison between the experimental and simulation results of *FM* of octorotor is shown in Figure 8.

**Figure 8.** Comparison between experimental and simulation results. 1.0 1.2 1.4 1.6 1.8 0.29

Spacing ratio

It can be seen that the simulation results are slightly higher than the experiment, with an error of about 6%, which may be the result of the introduced *C<sup>P</sup>* in experiments and the extra validation of *C<sup>P</sup>* for each simulation to reach the convergence.

#### *4.2. Simulation Results*

The velocity distribution between adjacent rotors with different rotor spacing is shown in Figure 9.

**Figure 9.** Velocity distribution. (**a**) *L*/*D* = 1.8; (**b**) *L*/*D* = 1.6; (**c**) *L*/*D* = 1.4; (**d**) *L*/*D* = 1.2; (**e**) *L*/*D* = 1.0.

According to Figure 9, it is not difficult to find that rotor wake began to attract each other as the rotor spacing was reduced. Since the wake of each rotor is affected by the wake of the two adjacent rotors, the wake of the rotor shows the phenomenon of inward contraction which may make the rotor unstable and affect the power consumption of the rotor. It also can be found that the vertical vortex below the side of the rotor slowly moves down as the rotor spacing decreases. In addition, the wake of 1.2 and 1.0 appears to be deeper than other spacing ratios, which may make the rotors more susceptible to ground effects.

The vortices distribution between adjacent rotors with different rotor spacing is shown in Figure 10.

**Figure 10.** Vortices distribution. (**a**) *L*/*D* = 1.8; (**b**) *L*/*D* = 1.4; (**c**) *L*/*D* = 1.0.

Figure 10 shows that the tip vortices and root vortices also attract each other and with a decrease in rotor spacing, the attraction between the vortices becomes more and more intense.

The pressure distribution around the rotor tips with different rotor spacing is shown in Figure 11.

**Figure 11.** Pressure distribution. (**a**) Pressure difference between upper and lower surface; (**b**) Pressure distribution of upper and lower surface.

As can be seen from Figure 11a, the pressure difference between the upper and lower surface at the blade tip is intended to decrease with a decrease in rotor spacing. It means a smaller thrust is obtained at smaller spacing. Figure 11b shows that the rotor spacing has a significant effect on the pressure on the upper surface of the rotor tip, but when the spacing is too small (such as *L*/*D* = 1.0), the pressure on the lower surface is also affected in this case.

The streamline distribution is shown in Figure 12.

The inward contraction of the rotor wake at 1.2 and 1.0 can be clearly observed from Figure 12d,e which indicates that the thrust and power variations are related to the location and the shape of the vortex. The Q-criterion for the octorotor with different rotor spacing is shown in Figure 13. As can be seen in Figure 13, the wake of the rotor shows the phenomenon of inward contraction at *L*/*D* = 1.0.

0.0 0.2 0.4 0.6 0.8 1.0

/ =1.0 / =1.4 / =1.8 -150 -120 -90 -60 -30 0 30 60

Pressure (Pa)

0.0 0.2 0.4 0.6 0.8 1.0

 lower surface-1.0 upper surface-1.0 lower surface-1.4 upper surface-1.4 lower surface-1.8 upper surface-1.8

/

Pressure difference (Pa)

**Figure 12.** Streamline for octorotor in hover. (**a**) *L*/*D* = 1.8; (**b**) *L*/*D* = 1.6; (**c**) *L*/*D* = 1.4; (**d**) *L*/*D* = 1.2; (**e**) *L*/*D* = 1.0.

**Figure 13.** Q-criterion for octorotor in hover. (**a**) *L*/*D* = 1.8; (**b**) *L*/*D* = 1.4; (**c**) *L*/*D* = 1.0.

#### **5. Conclusions**

To obtain the aerodynamic performance of a small octorotor with different rotor spacing, the small rotor test bench measurements were carried out with CFD simulations. The conclusions are as follows:


rotor. A decrease in rotor spacing will also affect the pressure difference between the upper and lower surfaces of the rotor tips. When the rotor spacing is too small, the pressure difference is sharply reduced, thus affecting the thrust generated. Future work should attempt to obtain an aerodynamic model that considers rotor interference to amend the conventional control theory with application in an octorotor with tilted rotors or wind disturbance.

**Author Contributions:** Y.L. carried out experiments, acquired the funding, and analyzed the results; Y.H. wrote the manuscript with assistance of Y.L.; H.W. performed a new set of simulations. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Nature Science Foundation of China (grant no. 51505087), the Fuzhou University Jinjiang Science and Education Park (no. 2019-JJFDKY-59) and the Fujian Provincial Industrial Robot Basic Component Technology Research and Development Center (2014H21010011).

**Acknowledgments:** The author thanks the Key Laboratory of Fluid Power and Intelligent Electro-Hydraulic Control (Fuzhou University), the Fujian Province University, and the Fuzhou University Jinjiang Science and Education Park for applying the experimental field.

**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*
