**1. Introduction**

Axial fans typically work in very turbulent flow conditions, e.g., because of their installations in pipelines, behind radiators, etc. This results in very unstable aerodynamic forces on the impeller blades, which in turn cause excessive sound radiation. Noise from flow machines consists of tonal noise, as a result of the interactions among the turbine blades and stationary housing components or guide vanes and broadband noise resulting from the acoustic signal generated by strong turbulent structures occurring in the flow. The most modern aeroacoustic computational methods enable increasingly reliable predictions of the generated noise. They usually require specific information about the transient flow field, obtained by simulation using computational fluid dynamics methods.

The most accurate of these methods, a direct numerical simulation (DNS), could solve the Navier–Stokes equation with no simplifications and could predict the unsteady flow and the associated acoustic field. Unfortunately, DNS is not feasible for a complex geometry, such as a fan, due to the enormous computational costs. Finding a non-stationary flow field with less effort requires modelling of at least part of the of the turbulent fluctuation [1]. Two different ways are currently used to reduce computational costs. The first is time averaging, which is known as the unsteady Reynolds-averaged Navier–Stokes (URANS) simulation (URANS), the second is a spatial filtering of the full Navier–Stokes equations, called large eddy simulation (LES).

In the case of URANS, the reduction of calculation costs is enormous but the cost is a large level of approximation. All random turbulent fluctuations are modeled, so only tonal sources of rotating machine sounds can be predicted. LES solves large turbulent structures, and only small eddies are modeled, but the computational costs are still high. Since this paper contains a very large number of numerical calculations that involve long-term calculations, the authors decided to use the URANS method in the simulations

**Citation:** Romik, D.; Czajka, I. Numerical Investigation of the Sensitivity of the Acoustic Power Level to Changes in Selected Design Parameters of an Axial Fan. *Energies* **2022**, *15*, 1357. https://doi.org/ 10.3390/en15041357

Academic Editors: Marcin Kami ´nski and Davide Astolfi

Received: 30 November 2021 Accepted: 7 February 2022 Published: 14 February 2022

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

**Copyright:** © 2022 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 (https:// creativecommons.org/licenses/by/ 4.0/).

carried out. The relevance of this decision is confirmed by Kissner et al. [2], who show that this way of modelling allows for satisfying accuracy to be achieved.

Solving aeroacoustics problems requires even more computational efforts. In addition to determining the sound sources, it is important to determine how the sound wave propagates. Assuming that the energy difference of the flow and the acoustic wave is sufficiently large, one can focus on the one-way coupling between the fluid flow and the acoustic signal generated. This line of reasoning led Lighthill to develop the aeroacoustic analogy named after him.

Lighthill [3–5] first made a formulation of the acoustic analogy for jet noise in 1952, which demonstrated that the flow mechanisms that were responsible for sound radiation could be expressed in terms of the quadrupole source. Curle [6] extended Lighthill's analogy to the fluid–structure interaction and implemented an extra acoustic source produced from the reaction force that was exerted on the fluids surrounding the body that did not move. FW-H [7,8] generalized Curle's analogy and extended the analogy to a moving structure. The FW-H equation splits up the aeroacoustic source into three different kinds of source: monopoles, dipoles, and quadrupoles.

The general theory of Lighthill introduced by FW-H, which takes into account the motion of a body as a potential source, has been used in many noise studies of rotating machines, i.e., turbines, fans, helicopters, etc. The FW-H analogy also takes surface sources into account, which makes it possible to determine noise from sources other than the quadrupole sources proposed by Lighthill. Schmitz and Yu [9] proved that, at a low Mach number, the volume integral makes no significant contribution to the noise generated by a hovering helicopter rotor. For a given range of speeds, the rotor is a surface of monopole and dipole sources and their contributions depend on factors, i.e., geometry, speed, and forces acting on the surface.

Brentner and Farassat [10], present a comprehensive review of the mathematical basis of the FW-H equation, comparing integral formulas and sufficiently powerful numerical methods applied to helicopter noise. They found that the contribution of the volume integral is small for subsonic flows, but gives a larger result for supersonic and transonic flows. In addition, they found that, by applying a permeable FW-H surface, instead of to the body surface, this would allow quadrupole sources to be included.

Konstantinov et al. [11] showed results of URANS, Delayed Detached Eddy Simulation (DDES), and LES flow and noise distributions in the test cabin segment. Compared to the FW-H mathematical model and the hybrid approach of solving the wave equation, including non-reflecting boundary conditions, a small influence on the sound pressure level from the imperfect boundary condition in the LES was shown.

Sundström et al. [12] used the LES method to investigate which acoustic sources predominated in low mass flux flows. They found that blade forces resulting from varying wall pressures are the main sources of generated noise at low mass flux flows. It turned out that sound sources coming from forces on blades (dipoles), were much larger than quadrupole sources, especially in subsonic flows. Comparing these two sound sources, one can see the relation *Wd* : *Wq* ∼ 1 : *<sup>M</sup>*<sup>2</sup> can be observed. Moreover, it can be concluded that monopole sources have a greater influence on the generated noise in sonic flows.

In a study by Al-Am et al. [13], the LES approach was used to numerically calculate the influence of selected parameters on the noise generated around a flat plate. Noise generated at the trailing edge and noise of turbulent nature was investigated. The flow character and geometry were chosen to correspond to the Amiet model. It is shown that the adopted model gives very good noise calculations in agreement with the analytical model and DNS calculations. Moreover, regarding the ACAT1 fan noise test [2], RANS-based analytical methods are commonly used to predict broadband fan noise. The accuracy of the aerodynamic noise results obtained in post-processing calculations depends not only on the choice of acoustic model itself, but to a large extent on the turbulent model adopted, which has a significant influence on the nature of the flow and which in turn affects the fan broadband noise. In continuation of this work, [14] focused on the importance of acoustic

models. Twelve different models were investigated and several different solvers were used to solve them. Both models—-based on acoustic analogies and those using direct methods—were compared. The methods used are distinguished by the turbulence models, the applied boundary conditions related to the propagation of the acoustic wave and the noise from the rotor blade phenomena. It turns out that at low frequencies, the differences in the generated noises are quite large, while at higher frequencies, the sound power level is within ±3 dB. Furthermore, it is proven that by increasing the rotational speed, the generated acoustic power is similar for different models.

Biedermann et al. in their study [15], provided detailed information on the broadband noise reduction possibilities of a low pressure axial fan with serrations on the leading edges. For the area of instability under partial load conditions, it is proposed that the dominant noise reduction mechanisms are dependent on aerodynamic effects related to the serrated geometry at the leading edge, which results in a reduction of the dominant low and medium frequency noise levels. It was also shown [16] that under the same operating conditions, the sound pressure levels at the two measurement points of the radial fan increases by approximately 5.8% and 2.8%, when the ambient pressure increases from 50 kPa to 100 kPa. As the ambient pressure increases, the fan sound pressure level shows an approximate logarithmic increase trend. It is worth mentioning that researchers [17] have attempted various techniques to identify the source of the sound, i.e., isocontours of the dilatation field, which revealed sources of acoustic scale disturbance, and may be the cause of the noise, a dynamic mode decomposition for the pressure upstream and downstream of the fan blade, which shows several strong fashions around the first three blade frequencies, and finally the acoustic analogy of FW-H, which showed a difference of about 5 dB between the blade tip and the lower parts of the blade in a specific frequency range. These results are consistent with the expectations that higher flow velocities would yield higher acoustic pressures.

The effect of the blade curvature on the generated noise was also investigated [18,19], in relation to the classic Amiet formulation. It has been proven that a curved blade causes a reduction in noise; this effect is particularly noticeable at the blade tip. Significant differences were also observed regarding noise generation at the leading and trailing edges of the blade, where the former is globally dominant but takes on values close to the latter at around 3.75 kHz and higher [20,21]. Similarly, other researchers confirm this relationship that the pressure fluctuation of the radial fan was smallest when the blade outlet angle was smaller, and it was also shown that a corresponding increase in the blade outlet angle reduced the amplitude of the pressure fluctuation in the blade pass frequency and its harmonics, which is conducive to reducing rotor noise [22]. In [23], it was shown that, on the upstream side on the blade walls, the sound pressure level was higher than on the downstream side. This was due to the separation of the boundary layer at the leading edge with increasing radial velocity near the ring, resulting in a low frequency noise. The leading edge therefore turned out to be the dominant dipole source generating tonal noise in contrast to the other rotor elements. The issues of stream separation at too low of a rotational speed of the rotor are also discussed [24], where under rotating stall conditions the fluctuation of the sound pressure amplitude becomes much greater than under other conditions, and the fluctuation of the sound pressure level is greater at a low frequency under stall conditions than under normal operating conditions.

The sensitivity analysis consists of examining how a given model depends on the parameters entering into it. By testing the sensitivity, you can determine which input variables have the greatest impact on a given output, and in this way, areas that require more attention may be identified [25]. Sensitivity analysis is a concept used in various engineering fields, i.e., acoustics, material science, environmental protection, etc. [26–28].

In the above works, different types of fans have been studied, both in terms of increasing their efficiency and reducing noise emissions, which proves the relevance of the issue. Researchers are still looking for a compromise to design these machines in such a way that they achieve the greatest efficiency with the least possible noise emissions. The authors of this paper, in previous studies [29–31], also discussed the influence of design parameters on the noise generated and the efficiency of different types of fans.

The research carried out in this work focuses on the sensitivity of the generated noise to changes in the design parameters of the fan rather than on the determination of the exact value of the acoustic power level.

### **2. Numerical Simulations**

#### *2.1. Research Object*

The numerical study was carried out on an axial fan (see Figure 1) with a diameter of 220 mm and a rotational speed of 3000 revolutions per minute, installed in a pipe of a circular cross-section with a diameter of 230 mm. The rotor has six blades on a hub with a diameter of 100 mm and a length of 200 mm at an angle of 20◦. For such an assumed rotational speed, the fan achieves a volume flow rate of 750 m3·h<sup>−</sup>1.

**Figure 1.** Geometrical model of the fan: (**a**) rotor, (**b**) pipeline.

The numerical model was a computational grid consisting of 2,120,050 cells, in which the Navier–Stokes equations were solved using the finite volume method. The main fan parameters are shown in Table 1.



In the model, two characteristic zones can be identified: a rotating zone in which the rotor is located and a stationary zone representing a curved pipeline. The geometric model does not include the fan fixing elements inside the pipeline. To simplify the model, steering systems and airflow straightening elements were also omitted. Marked in Figure 2, distances *Lp* and *Lt* are equal to each other and are 70.16 mm. Inlet straight sections *Li* and outlet *Lo* of the pipeline are equal and amount to 100 mm.

**Figure 2.** Diagram of built-in fan.

The geometric model of the rotor was made in SolidWorks and then adapted to the ICEM CFD environment. The ICEM CFD program was used to create a orthogonal computational mesh based on a multi-block structure grids topology. In the areas close to the leading and trailing edges, vortices were created due to the separation of the boundary layer, which is why the computational grid was appropriately compacted at these points in order to solve the flow with greater precision. In addition, a wall layer was modeled to provide a parameter *y*<sup>+</sup> < 5. In order to analyze the validity of the model performed, an analysis of the independence of the results from the calculation grid was performed. The mesh used (see Figure 3) allowed for results of sufficient accuracy to be obtained in a reasonable time.

**Figure 3.** Fan computational grid: (**a**) rotor mesh, (**b**) pipeline mesh.

#### *2.2. Mathematical Model*

The first and very important step in the simulation is to locate the sources that can be used to calculate the generated noise. This can be achieved using appropriate computational fluid dynamics techniques. In this field, the basic equations are based on the Navier–Stokes equations derived from the conservation of mass (continuity equation), conservation of momentum, and conservation of energy.

The continuity equation can be written for fluids as:

$$
\rho \frac{\partial u\_i}{\partial x\_i} = 0 \tag{1}
$$

where *ρ* is the density of the fluid, *t* is the time, components of the velocity vector **u** in the coordinate system, and *xi* coordinates in the Cartesian system. The momentum equation is written as:

$$
\rho \frac{\partial u\_i}{\partial t} + \rho \frac{\partial u\_i u\_j}{\partial x\_j} = -\frac{\partial p}{\partial x\_i} + \frac{\partial v\_{ij}}{\partial x\_j} \tag{2}
$$

where *p* is the pressure and *τij* are viscous stresses. Equation (2) is derived by applying Newton's second law of dynamics, which relates the forces acting on a fluid volume to its acceleration.

Turbulence model *k* − *ω* is one of the most popular models, which shows the phenomena of turbulent flow. It belongs to a family of models in which all turbulence effects are modeled. This is a two equation model. This means that the transport equations are solved to include phenomena, such as convection and turbulent energy dissipation. The variables considered in the equation are the turbulent kinetic energy *k*, representing the turbulence energy, and the specific turbulent dispersion coefficient *omega*, denoting the dispersion rate of the turbulence kinetic energy. Variable *ω* is also known as the turbulence scale. The standard *k* − *omega* model works well for low Reynolds number flows where the boundary layer is appropriately sized and the viscous sublayer is well separated. Thus, the standard model *k* − *ω* is best suited for modelling the boundary layer. Other advantages include excellent performance in complex near-wall flows with adverse pressure gradients and separation, e.g., in rotating machinery. The model also predicts excessive and early vortex separations.

Model *k* − *ω* SST is a model that offers the strengths of the *k* − proposed by Launder and Spalding [32] and model *k* − *ω* proposed by Wilcox [33], and provides an additional component to limit the overproduction of turbulent kinetic energy in areas of high pressure gradients (stagnation points, areas of separation vortex near wall layer). Menter [34] examining models *k* − and *k* − *ω*, and observed that the first handles turbulence well in free and shear layers and shows negligible sensitivity to inlet boundary conditions for quantities describing turbulent flow. This is a desirable quality because these quantities are often not exactly known in practical calculations. However, the *k* − *ω* model better models turbulent flow in the boundary layer but is more sensitive in free flow.

The sound pressure level (SPL) was determined using the FW-H analogy. This model is based on the Lighthill analogy and allows noise to be determined by equivalent acoustic sources. Ansys Fluent uses these equations to determine the sound pressure at a given distance from a sound source by an integral over the surface containing those sources. The FW-H equation is a non-homogeneous wave equation [7,35], which can be derived by combining the continuity and Navier–Stokes equations. It can be written as

$$\begin{split} \frac{1}{\alpha\_0^2} \frac{\partial^2 p'}{\partial t^2} - \nabla^2 p' &= \frac{\partial^2}{\partial x\_i \partial x\_j} \{ T\_{ij} H(f) \} \\ &- \frac{\partial}{\partial x\_i} \left\{ \left| P\_{ij} n\_j + \rho u\_i (u\_n - v\_n) \right| \delta(f) \right\} \\ &+ \frac{\partial}{\partial t} \left\{ \left| \rho\_0 v\_n + \rho (u\_n - v\_n) \right| \delta(f) \right\} \end{split} \tag{3}$$

where *ui*—air velocity in the direction of *xi*, *vi*—surface velocity in the direction of *xi*, *un*—air velocity normal to the surface *f* = 0, *δ*(*f*)—Dirac delta, *vn*—velocity of the surface normal to the surface, *H*(*f*)—Heaviside function, *p* —sound pressure in the far field (*p* − *p*0), *ni*—normal vector pointing to the external area (*f* > 0), *a*0—speed of sound in the far field, *Pij*—compressive stress tensor, *Tij*—Lighthill stress tensor.

To solve Equation (3), the Green's function must be used to the open area. The complete solution involves the calculation of surface and volume integrals, the first representing monopole, dipole, and partially quadrupole acoustic sources, and the second representing quadrupole sources in the area outside of the source surface. The volume integral becomes negligible when the Mach number value of the flow is small and the source area covers the source area. In Ansys Fluent, choosing a source on a solid surface-like rotor, the volume integrals are neglected, then the equation takes the following form:

$$p' = (\vec{\mathbf{x}}, t) = p'\_T(\vec{\mathbf{x}}, t) p'\_L(\vec{\mathbf{x}}, t) \tag{4}$$

$$\begin{split} 4\pi p\_T'(\vec{x}, t) &= \int\_{f=0} \left[ \frac{\rho\_0 (\mathcal{U}\_n + \mathcal{U}\_h)}{r (1 - M\_r)^2} \right] dS \\ &+ \int\_{f=0} \left[ \frac{\rho\_0 \mathcal{U}\_h \left\{ r \mathcal{M}\_r + a\_0 \left( \mathcal{M}\_r - \mathcal{M}^2 \right) \right\}}{r^2 (1 - M\_r)^3} \right] dS \end{split} \tag{5}$$
 
$$\begin{split} 4\pi p\_L'(\vec{x}, t) &= \frac{1}{a\_0} \int\_{f=0} \left[ \frac{L\_r}{r (1 - M\_r)^2} \right] dS \\ &+ \int\_{f=0} \left[ \frac{L\_r - L\_M}{r^2 (1 - M\_r)^2} \right] dS \\ &+ \frac{1}{a\_0} \int\_{f=0} \left[ \frac{L\_r \left\{ r \mathcal{M}\_r + a\_0 (M\_r - M^2) \right\}}{r^2 (1 - M\_r)^3} \right] dS \end{split} \tag{6}$$

where

$$\mathcal{U}\_i = v\_i + \frac{\rho}{\rho\_0} (u\_i - v\_i) \tag{7}$$

$$L\_i = P\_{i\bar{j}}\hat{n}\_{\bar{j}} + \rho u\_i (u\_{\bar{n}} - v\_{\bar{n}}) \tag{8}$$

The contribution of quadrupole terms (volume integrals) in the FW-H analogy is proportional to the square of the Mach number (*M*2). In the analyzed system, the Mach number reaches values below 0.1, which means that the volume integrals can be omitted. Considering the time *t* and a distance to the observer *r*, the integral equation takes into account the delay due to the distance from the source to the receiver, according to the following formula:

$$
\pi = t - \frac{r}{a\_0} \tag{9}
$$

$$\begin{array}{ll} M\_r = M\_i r\_i & \dot{M}\_r = \frac{\partial M\_i}{\partial \tau} r\_i \\ Q\_{\text{il}} = Q\_i \dot{\boldsymbol{n}}\_i & \dot{Q}\_{\text{\textell}} = \frac{\partial Q\_i}{\partial \tau} \dot{\boldsymbol{n}}\_i & Q\_{\text{\textell}} = Q\_i \frac{\partial \boldsymbol{n}\_i}{\partial \tau} \\ L\_i = L\_i \dot{\boldsymbol{n}}\_i & L\_r = \frac{\partial L\_i}{\partial \tau} \dot{\boldsymbol{r}}\_i & L\_r = L\_i \boldsymbol{\mathfrak{f}}\_i & L\_M = L\_i M\_i \end{array} \tag{10}$$

where *n*,*r*—unit vectors of radiation and normal to the wall, *M*—Mach number of the surface source velocity component along the direction of the radiation vector.
