3.1.4. *RNG k*−ε Turbulence Model

The *RNG k*−ε model is derived from the mathematical method of renormalization group in combination with the transient Navier-Stokes equations (N-S equations). This model is similar to the standard model and its analytical characteristics directly evolved from the standard model. The main difference between the RNG model and the standard one is due to the consideration of the turbulent vortex by the RNG model. This consideration enhances the calculation precision on the vortex. Moreover, the turbulent Prandtl number also provides a comprehensive analytical equation. A condition is also included into the turbulent diffusion equation in order to improve the precision of the standard model. The flow field can be presented in a more precise way. The equation of the RNG model is described as follows [18].

Equation of turbulence kinetic energy (*k*)

$$\frac{\partial}{\partial t}(\rho k) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho k u\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \Big(a\_k \mu\_{eff} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_j}\Big) + \mathbf{C}\_k + \mathbf{G}\_b - \rho \varepsilon - \mathbf{Y}\_M \tag{11}$$

Equation of diffusivity (ε)

$$\frac{\partial}{\partial t}(\rho \varepsilon) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho \varepsilon u\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \Big( a\_k \mu\_{eff} \frac{\partial \varepsilon}{\partial \mathbf{x}\_j} \Big) + \mathbb{C}\_{1\varepsilon} \frac{\varepsilon}{k} + (\mathcal{G}\_k + \mathcal{C}\_{3\varepsilon} \mathcal{G}\_b) - \mathcal{C}\_{2\varepsilon} \rho \frac{\varepsilon^2}{k} - \mathcal{R}\_{\varepsilon} \tag{12}$$

where α*<sup>K</sup>* and αε are the turbulence Prandtl number of the turbulence kinetic energy and the turbulence diffusivity, μ*eff* is the coefficient of equivalent turbulence viscosity, *R*<sup>ε</sup> is the parameter of the modified turbulence viscosity, the constants are *C*1<sup>ε</sup> = 1.42, *C*2<sup>ε</sup> = 1.68 respectively.

The main difference between the RNG model and the standard model is described as follows. The RNG model is used to build up a new equation based on the condition of low Reynolds number. The equation is as follows.

$$d\left(\frac{\rho^{2k}}{\sqrt{\varepsilon\mu}}\right) = 1.72 \frac{\stackrel{\frown}{\upsilon}}{\stackrel{\frown}{\upsilon}^3 - 1 + \text{C}\_V} \tag{13}$$

where *CV* <sup>≈</sup> 100, <sup>∩</sup> *<sup>v</sup>* <sup>=</sup> <sup>μ</sup>*eff* <sup>μ</sup> .

The equation depicts how the Reynolds number affects the coefficient of equivalent turbulence viscosity so that a model could perform better at a low Reynolds number. At conditions with a higher Reynolds number, the turbulence velocity equation of the standard model is still used, except that the *C*μ parameter is set as 0.0845 according to the RNG theoretical calculation.

Moreover, since the turbulence in a uniform flow is also affected by the vortex, the turbulence viscosity is also modified to compensate for this influence. The equation is as follows.

$$
\mu\_t = \mu\_{t0} f(\alpha\_s, \Omega, \frac{k}{\varepsilon}) \tag{14}
$$

Here μ*t*<sup>0</sup> is the quantity that is not modified from the original equation of the turbulence viscosity coefficient. Ω is a characteristic parameter that is used by the FLUENT software. α*<sup>s</sup>* is a vortex constant and its value is determined form the vortex intensity of the flow condition. For a moderate vortex flow state, α*<sup>s</sup>* is set as 0.05. For a stronger vortex flow, a larger value can be used.

For the turbulence Prandtl number, the RNG theory supplies a comprehensive analytical equation that can be used to calculate α*<sup>k</sup>* and αε as follows.

$$
\left|\frac{\alpha - 1.3929}{\alpha\_0 - 1.3929}\right|^{0.6321} \left|\frac{\alpha + 2.3929}{\alpha + 2.3929}\right|^{0.3679} = \frac{\mu\_{mol}}{\mu\_{eff}}\tag{15}
$$

where α<sup>0</sup> = 1.0 and α*<sup>k</sup>* = αε ≈ 1.393 for a larger Reynolds number.

Finally, the conditional parameter of *R*<sup>ε</sup> is also included into the diffusivity equation. This parameter leads to the main difference from the standard model and the equation is as follows.

$$R\_{\varepsilon} = \frac{\mathbb{C}\_{\mu} \rho \eta^{3} (1 - \eta / \eta\_{0})}{1 + \beta \eta^{3}} \frac{\varepsilon^{2}}{k} \tag{16}$$

where η = *Sk* <sup>ε</sup> , η<sup>0</sup> = 4.38, β = 0.0 = 0.012.

Since the RNG model provides a comprehensive definition and correction to several parameters, the RNG model can react to the flow field with immediate changes and curved streamlines. This is also the reason why the RNG can present better performance in this type of flow field.

#### *3.2. Performance Testing Equipment for Fans*

From the aspect of performance measurements in this experiment, the detailed configuration of the measuring equipment and apparatus and the corresponding operations were described separately as shown in Figure 4. The main device of the performance testing equipment for the fans used adopted the outlet-chamber wind tunnel according to the AMCA 210-99 standard. This type of wind tunnel is

composed of a main body, flow setting means, multiple nuzzles, and a flow-rate regulating device [14]. Its major function is to simulate various types of the air-flow condition downstream of the fan and supply a good and stable flow field for measurement in order to obtain the complete performance curves [26].


11 Adjusting device for fiber-optic tachometer 22 Access hole

**Figure 4.** Configuration of the fan performance testing system and specifications of the major instruments.

#### 3.2.1. Calculation of Flow Rates

According to the calculation by the standard equations and flow-rate measurements by the National Laboratory, the errors can be obtained through a comparison with the measurement readings. The standard equations are as follows [27].

The pressure difference between the nozzle outlet and inlet *PL*<sup>5</sup> and *PL*<sup>6</sup> can be obtained. The flow rates on the cross-sections of the nozzles can be determined with varying nozzle coefficients as shown in Figure 5. If there is a need to calculate the outlet flow rate of the fan under test, then the effect of density variations must be considered, the measurement of which is as follows [28].

**Figure 5.** Definition of various measurement planes.

The equation for the calculation of flow rates in a test chamber with multiple nuzzles is

$$Q\_{\mathfrak{F}} = 265.7Y\sqrt{\Delta P/\rho\_{\mathfrak{F}}} \sum\_{n} \left( \mathbb{C}\_{n} A\_{\theta n} \right) \tag{17}$$

where *Q*<sup>5</sup> is the total flow rate measured by a bank of nozzles, CMM; Δ*P* is the pressure difference across the nozzles, mm-Aq; ρ<sup>5</sup> is the air density upstream of the nozzles, kg/m3; is the expansion factor; *Cn* is the discharge coefficient of the *n*th nozzle; and *A*6*<sup>n</sup>* is the cross-sectional area of the *n*th nozzle's throat, m2.

#### 3.2.2. Calculation of Air Pressures

Typical pressure readings can be directly measured by instruments, but requires understanding the definition of a fan's static pressure (Δ*Ps*) and total pressure (Δ*Pt*). The static pressure, defined as the difference between the fan's static pressure at typical pressure readings, can be directly measured by instruments, but understanding (*Ps*<sup>2</sup> ) and the static pressure at inlet (*Ps*<sup>1</sup> ) is required. The total pressure is the difference between the fan's total pressure at outlet (*Pt*<sup>2</sup> ) and the total pressure at inlet (*Pt*<sup>1</sup> ). The equations for measurement and calculation are explained respectively as follows.

Since the outlet and inlet planes of the fan under test are *PL*<sup>2</sup> and *PL*1, respectively, they can be defined as follows [29]:

$$P\_s = P\_t - P\_v \tag{18}$$

$$P\_t = P\_{t\_2} - P\_{t\_1} \tag{19}$$

where *Ps* is the static pressure of the fan under test; *Pt* is the total pressure of the fan under test; *P*<sup>ν</sup> is the dynamic pressure of the fan under test; *Pt*<sup>2</sup> is the total pressure at the fan's outlet (or plane *PL*2); and *Pt*<sup>1</sup> is the total pressure at the fan's inlet (or plane *PL*1).

Since in this experiment there was no duct at the inlet of the fan under test, therefore *Pt*<sup>1</sup> = 0. On the other hand, the measured static pressure at the outlet is the same as the static pressures measured at the measuring plane *PL*7. Therefore, *Ps*<sup>2</sup> = *Ps*<sup>7</sup> .

$$P\_{t\_2} = P\_{s\gamma} + P\_v \tag{20}$$

$$P\_{\ $} = P\_{\$ \gamma} \tag{21}$$

It can be concluded from the above equation that the static pressure of the fan under test happened to be equal to the static pressure obtained at the outlet test chamber *Pt*<sup>7</sup> . The Type A method (a test method with no duct at either the outlet or the inlet) that was carried out at the outlet test box is a special case of testing. When carrying out different types of tests or different equipment, the equation of *Symmetry* **2019**, *11*, 959

the static pressure of the fan under test is thus more complicated. The calculation of dynamic pressure is as follows [30]:

$$P\_{v\_2} = \frac{\rho\_2 V\_2^2}{19.6} \tag{22}$$

where *P*ν<sup>2</sup> is the outlet dynamic pressure of the fan under test, mm-Aq; *V*<sup>2</sup> is the outlet air velocity of the fan under test, m/s; ρ<sup>2</sup> is the outlet air density of the fan under test, kg/m3; and *<sup>V</sup>*<sup>2</sup> = *<sup>Q</sup>*<sup>2</sup> <sup>60</sup>*A*<sup>2</sup> <sup>=</sup> *<sup>Q</sup>* <sup>60</sup>*A*<sup>2</sup> · <sup>ρ</sup> ρ2 = *<sup>Q</sup>* <sup>50</sup>ρ2*A*<sup>2</sup> where *Q*<sup>2</sup> is the outlet flow rate of the fan under test, CMM; *Q* is the standard flow rate of the fan under test, CMM; *A*<sup>2</sup> is the outlet cross-sectional area of the fan under test, m2; ρ is the density of air at STP (1.2 kg/m3); and *Pt* = *Ps* + *Pv* = *Ps* + *Pv*<sup>2</sup> .

$$P\_t = P\_s + \frac{\rho\_2 V\_2^2}{19.6} \tag{23}$$

#### 3.2.3. Fan Performance Power and Efficiency

The calculation of power can be obtained from torque and rotation speed. By measuring the torque of a fan by a torque gauge and measuring the rotation speed by a fiber-optic tachometer, the fan input power (*W*) can be obtained. The fan efficiency can also be obtained from the air pressure and the flow rate. It can be estimated by the equations as follows:

$$\mathcal{W} = \frac{2\pi \times T \times n}{33,000 \times 12} \tag{24}$$

$$
\eta\_s = \frac{P\_s \cdot Q}{4500 \text{W}} \tag{25}
$$

$$
\eta\_t = \frac{P\_t \cdot Q}{4500\,\text{W}}.\tag{26}
$$

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

The flow passage of the numerical model of the target axial fan is shown in Figure 6. Within the entire range of the flow rate, the model was composed of not only one flow passage. It is known that the flow field behavior follows the continuity equation and the momentum equation. The axial fan is composed of the inlet cone, impeller, and the fan housing. The resulting mesh structure contains structured meshes as the majority and unstructured meshes as the minority. After further mesh refinement and coupling, the final number of the mesh contained 750,000 cells as shown in Figure 7. The mesh structure contained a rotating mesh system between the inlet and the outlet. The outlet boundary had a uniform pressure of 1 atm. The rotation speed of the rotating mesh was 2000 RPM. The pressure and velocity coupling wind-facing difference method was selected as the simple algorithm. The maximum residual was defined as <10−<sup>3</sup> for convergence.

Since this is a problem for rotating machinery, the commercial CFD software FLUENT6.3 was used for the simulation. The selected region was configured to rotate against an axis so the momentum equation could be rectified automatically. The source terms were automatically added into the relevant equations for calculation. The configuration of the boundary conditions needs to consider the operating condition of a real object, i.e., to comply with the physical phenomenon. Otherwise, the accuracy of the calculated result might be affected. The boundary conditions of this simulation include the inlet boundary condition, outlet boundary condition, and wall boundary condition. Their descriptions are listed in Table 4 as follows as shown in Figure 8 [15,28,31].

**Figure 6.** Numerical model of the target axial fan.

**Figure 7.** Mesh structure along the cross-section of the numerical model of the axial fan.

**Table 4.** Boundary conditions of the fan model for simulation.


**Figure 8.** Boundary conditions.

### **5. Verification of the Case Study and Numerical Analysis**

#### *5.1. Verification between Numerical Simulation and Experiment Testing*

After the completion of the numerical simulation, we compared the results to the experiment in order to determine whether the trend corresponded to each other. Furthermore, the accuracy of the numerical simulation in this study was also verified as shown in Figure 9. First, to verify the accuracy of the numerical simulation, it is known from the comparison of the numerical results to the experimental ones in Table 5 that the numerical results were on average 3% smaller than the experimental ones. After exploring the reason, it was found that the deviations occurred most from the difference between the configuration of the real test environment and that of the simulated environment. This is due to the fact that when testing a real fan, the air impedance varies according to different operating points. It is also known from the experimental results that the average fan speed during the real tests was 3% larger than the design fan speed, so that it could be closer to the real test result [34].

**Table 5.** Comparison of the numerical simulation result and the experimental result.


By comparing the numerical simulation results to the experimental results, the non-dimensional air velocity at the inlet, as shown in Figure 10, served as a good reference for analyzing the flow field distribution of these six opening patterns. From the average value at the inlet and the pulse velocity distribution, the correct air velocity can be determined from Figure 10. The results indicated that the deviation of the non-dimensional velocity at the inlet 4H was the same for the six opening patterns on the boundary. Therefore, the mesh structure was valid for the simulation in this study.

**Figure 9.** *Cont.*

**Figure 9.** Resulting plots of CFD simulation for the numerical fan model. (**a**) Pressure distribution on top and bottom covers; (**b**) pressure distribution on the fan; (**c**) distribution of stream lines; (**d**) distribution for flow field; and (**e**) program convergence plot.

Figure 11 indicates the non-dimensional velocity V along the vertical direction from A to I. Since the air flow rate of a fan is affected by the gaps at the outlet, the relative location of the intake channel is different according to its air velocity. The velocity component of the six curves on 4H in Figure 11 was consistent with those in Figure 10. The velocity component V at 10 mm in front of the fan was almost uniform. The variation in the air velocity was due to the interference of the direction of the rear gaps. At a distance of −1H, the velocity started to vary and the change in the air velocity was in the same direction. A similar trend could be observed at the eight vertical lines at the rear end of −3H. The V component was larger at both the left and right sides and was negative at the centerline of the fan. It is known that the magnitude of the velocity component at −5H is relatively smaller than the left and right sides. The results also indicated that the velocity component of the Idea-D magnitude remained at the highest air velocity from the measurement line between E and I. Under this condition, the boundary layer effect of the velocity component could be clearly observed. The magnitude of the velocity component was gradually increased and the Idea-D can be viewed as the optimal factor of the opening pattern design.

**Figure 10.** Non-dimensional velocity V/Uref at nine planes along the vertical centerline of the model for the numerical simulation.

ŗİŖųŦŧ

**Figure 11.** Non-dimensional velocity V/Uref at nine planes along the vertical centerline of the model for the numerical simulation for these six opening patterns.

#### *5.2. Design Cases and Comparison between Simulation Results*

Figure 12 indicates the qualitative velocity distribution on the three-dimensional plane at location F of the outlet opening. The purpose of this figure is to understand the velocity pattern on these six planar regions. As the results clearly indicated a cone-shaped distribution, the velocity distribution at the outlet region was the most apparent. Moreover, Figure 12 also indicates that the out ring of the distribution of Idea-C revealed a low air velocity. In contrast, the air velocity at the peak was the highest. The Idea-D had the largest number of peaks and this indicates that this region was the region with the highest air velocity and the distribution had a more focused region in the central region. For the velocity distribution of Idea-B, E, and F, since they were affected by the non-uniform gaps, the air flow velocity was reduced. When the fan was rotating, the air layers also rotated. The air at the central region was affected by the viscosity and molecular attraction, so the rotating effect was transferred to the air on the out ring of the flow field.

**Figure 12.** Three-dimensional distribution of the air velocity at the outlet.

The results of the simulation for the six opening patterns are shown in Table 6. The effect of different opening patterns can be evaluated by a comparison of the numerical results. The variation in the flow field of six different opening patterns and the change in the performance were reviewed from the inlet and outlet of each streamline. Table 6 reveals that the high velocity regions were on the left and right sides of the flow field. This phenomenon also directly affected the outflow of the air between the blades. This is the same phenomenon that reduced the range of higher air velocity on the downstream of the fan. For Idea-A, the expansion of the opening pattern provided some improvement and became more and more apparent. The region of low air velocity extended inward so that the velocity along the inner diameter increased gradually. For Idea-B and E, the region of the highest air velocity at the outlet reduced gradually. The region of low air velocity also expanded along the downstream direction. The downstream flow field of the Idea-D honeycomb type opening still remained smooth, but not affected by the opening gap. The flow field revealed high velocity at the middle portion and downstream portion. On the overall distribution, the outlet of the impeller still remained at the high air velocity.

**Table 6.** Comparison of the numerical simulation results for six different cases.

Moreover, the flow field of those six opening patterns indicated a trend of moving downstream toward the outlet. A non-ideal velocity distribution at the outlet could severely affect the total air flow rate. It is known from Table 6 that Idea-B presented a very non-uniform velocity distribution. Most of the fluid aggregated to the left and right sides of the circular region. Reverse airflow and low air velocity can be found at the central region. The total air flow rate of the fan was not greatly affected by this phenomenon, and therefore this portion could be excluded.

The comparison of the resulting flow rate of the six opening patterns is shown in Figure 13. At the outlet of Idea-D with the honeycomb shape, the flow rate obtained from the numerical simulation was 18CFM. Similarly, the resulting flow rate of Idea-A was 16CFM. Figure 13 also reveals that the air flow rate of Idea-D was 2CFM larger than Idea-A. For the static pressure, the A4 model had the highest static pressure of 1.87 mm-H2O. The result static pressure of Idea-A was 1.74 mm-H2O. From the aspect of overall air flow rate, the simulation results indicated that Idea-D had the largest air flow rate.

**Figure 13.** Comparison of air flow rate and static pressure of new fan designs by simulation.

#### **6. Conclusions**

In this study, the aerodynamic performance of axial fans was investigated to determine the effect of different opening patterns. As the air flow rate decreases with increasing pressure, the fan's aerodynamic performance is affected if the density of the opening pattern is increased, so the boundary layer is affected. The cross-sectional flow field of an axial fan was reviewed and investigated. The flow field could be improved by making the streamlines smoother. Some reverse flows were also observed at the outlet of the fan. Conclusions of the investigation of the flow field along the radial direction and at the outlet are as follows.


There are many design parameters when desigrning the opening patterns of a fan outlet. The analysis in this study investigated six factors in order to provide future fan designs with a good reference for the optimization of the opening pattern. The results indicate that the honeycomb-shaped opening pattern can be beneficial to the performance enhancement of an axial fan. The parameters of the inlet design can also be included for further investigation on design optimization.

**Author Contributions:** The author contributed to the paper. H.-H.L. Collects and organizes data and acts as the corresponding author, J.-H.C. and the authors propose methods.

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

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