Reducing the noise produced by the interaction between a turbulent subsonic flow and a solid’s surface can be difficult, especially in industrial configurations where the noise level is a quality and selection criterion. When considering subsonic turbomachines, this noise can be reduced if the geometry of the profile is appropriate. However, there is no general solution, and each fan configuration needs a specific design. Therefore, most of the time acoustic engineers use a case by case approach, either experimentally or numerically   . In this context, identifying and characterizing the region responsible for the noise interaction becomes a challenge. Since the acoustic source is connected to the flow some studies for example separated the turbulent field into acoustically radiating and non-radiating components using Navier-Stokes equations   . Other approaches are based this time on the decomposition of the sound radiated itself in order to identify the most radiating zones as the case in our investigation.
This paper proposes a method of identifying and analyzing these zones based on the combination of Ffowcs Williams and Hawkings’ acoustic analogy (FW&H)  and either proper orthogonal decomposition (POD) or singular value decomposition (SVD)  . The SVD approach used in this paper is the application of the POD to a matrix formed by the product of the correlation matrix obtained from the elementary acoustic pressure to a receiver and it transposed. This study will develop the POD approach in more detail.
The pressure fluctuations obtained from a CFD approach are used as input data for the FW&H acoustic analogy  . For that, the large eddy simulation (LES) method   proposed by OpenFOAM  is used. OpenFOAM’s popularity has been growing rapidly for several applications, with an increasing number of users over the last few decades. Numerous investigations have been published in the CFD field. For example, de Villiers  investigated the influence of subgrid scale models in the method implemented in OpenFOAM. This investigation was carried out on several flow configurations around a rigid solid. The noise generated by the interaction between an airflow and the rearview mirror of a car was also studied. The results are in agreement with the experimental results. Furthermore, Dimtry et al.  modeled the flow around a cylinder with Smagorinsky’s subgrid models and the one-equation turbulent kinetic energy model implemented by de Villiers. The results are in agreement with experimental data   . Several other studies in the literature have demonstrated OpenFOAM’s effectiveness in modeling complex configurations   . In acoustics, in addition to de Villiers’ work, Olivier’s work  addressed the noise of the trailing edge of NACA 0012. It used data from the OpenFOAM CFD to supply the acoustic source of the Curle analogy. These acoustic results are in agreement with Herr et al.’s results  . The subgrid model of one of de Villiers’s turbulent kinetic energy equations will be used for this study for the CFD.
POD was first introduced by Lumley in 1967  in fluid mechanics as an objective method to identify and extract coherent turbulent flow structures. Building on this idea, Bonnet et al.  used POD in conjunction with a stochastic estimation method to identify the key characteristics of a turbulent field on the basis of a reduced number of measurements. POD was also employed either for optimal flow control   or reduction model development   . In the latter application, the Navier-Stokes equations are projected on the eigenmodes using the Galerkin projection. The POD method was used in many aeroacoustics investigations, mostly to describe the noise radiated by turbulent flows (e.g. jet noise, wake noise)   . The optimal decomposition basis velocity field is then used to estimate acoustic noise. Druault et al.   employed this method in the far field region to separate the acoustic contribution of the most energetic structures from the residues caused by the clean flow noise (i.e. the noise from velocity fluctuations). They noticed that 99.8% of the acoustic energy comes from the dominant modes and 0.2% from acoustic residues. Hekmati and Ricot  reached the same conclusion by applying the Druault et al.’s method   for wake noise generated by the blades of an axial fan. Furthermore, Gleggs and Devenport  used POD to describe the input to the turbulence using a set of statistically independent modes of velocity. The wall pressure is linked to the input turbulence for the evaluation of the acoustic interaction sound pressure.
SVD is generally used to search for the propagation operators such as the Green function   , or to solve inverse problems in acoustics     . To the authors’ knowledge, SVD has not yet been used to locate radiating areas on a moving surface due to its interaction with flow.
In the light of previous works, this study applies the POD and SVD methods to the problem of aeroacoustic noise generated by the interaction between a stationary blade and a turbulent flow in a channel. The objective is to understand the link between the decomposition modes of POD or SVD and the noisiest zones of the blade’s surface. The methodology relies on a three-step methodology: 1) The internal flows of the centrifugal fan are modeled using LES method   . OpenFOAM extend 3.2  software is used. The objective is to estimate the wall pressure fluctuation on the blade. 2) The previous wall pressure fluctuation of the blade is used to estimate the loading noise from the FW&H analogy   . 3) Finally, POD and SVD are used to extract the most important acoustic modes and visualize them in order to identify the zones that radiate the most on the blade’s surface. This paper is organized as follows. First, the estimation of the acoustic field based on the FW&H analogy is introduced. Then both the POD and SVD methods are developed. Finally, an application is demonstrated using a stationary blade.
2.1. Estimation of the Acoustic Field (FW&H)
The acoustic approach is based on the FW&H analogy  and is used as the reference approach. The subsonic regime is considered in this paper. The finite thickness of the blades is neglected; only loading noise ( ) generated by the fluctuating wall pressure on the blade is considered in the following example.
To alleviate the computation time problem, Formulation 1A proposed by Farassat  can be used. In this approach, the receiver time derivative in Formulation 1  is transformed into a retarded time derivative. This has also the great advantage of permuting the time derivative and the integral. Moreover, no derivation of the integral is needed. For the case in which the source and the receiver are stationary and the propagation medium is at rest, the Farassat Formulation 1A  in the far field becomes:
where is the source surface, is the distance between the source position on the blade’s surface and the receiver position , c is the speed of sound of the acoustic medium at rest, with is the unit normal vector to the source’s surface, and p is the wall pressure fluctuation of the blade’s surface obtained by CFD calculation. is the Mach vector number and and indicates that all the integrands should be evaluated at the retarded time , with receiver time t.
As mentioned previously, Formulation 1, proposed by Farassat, has the main advantage of avoiding the spatial derivatives; however, the receiver time derivative is maintained on and . Implementing this operation is complex, and the computation time increases. To evaluate Equation (1), two computational approaches are available in the literature   using either retarded time or advanced time. The latter is usually chosen when the aerodynamic data comes from CFD computations, as in this study. The received acoustic pressure must then be determined while an irregular receiving time discretization appears despite the regular emission time. An interpolation is thus necessary to obtain the received sound pressure at regular time.
The advanced time approach and the Lagrange interpolation are used in this paper. The calculation of the sound pressure at receiver in the far field and in a free medium allows for the calculation of root mean square sound pressure. When considering the loading noise (Equation (1)), the mean square acoustic pressure reads:
where is the temporal average over the time period . Thus, the radiated sound power estimated from far field microphones (i.e., receivers) located on a spherical surface encompassing the source can be written:
where is the receiving surface, is the elementary surface associated with receiver , and is the density of the surrounding fluid medium.
2.2. Proper Orthogonal Decomposition (POD) Approach
Generally, POD used in aeroacoustics is not based on acoustic analogies. The approach developed in this paper is a combination of the dipole term (loading noise) of the FW&H analogy and the POD theory. If one considers a dipole located at , the sound pressure received at point at time t according to Equation (1) can be written:
where is the i-th elementary surface of the source located at . Here, this elementary surface comes from the discretization of the source’s surface for the LES calculation. Since the sound pressure received at one point corresponds to the contribution of all the dipoles located on the wall surface, two matrices and are defined by
where represents the emission time for the dipole source located at at time step . Each column vector of the matrix is the sound contribution of all dipoles at a given reception time, while each row vector represents the sound pressure of one single dipole source along the receiving time. Matrix in Equation (5) is the correlation matrix of the sources for receiver . It is symmetric, real, positive definite, and spatial. Its eigenvalues (modes) are thus real, positive, and space-dependent. When developing a modal basis for matrix , one defines and , the diagonal matrix of eigenvalues and the matrix of eigenvectors at receiver respectively. Each column is the eigenvector associated with the eigenvalue at receiver . Then for every receiver , the problem to be solved is the following eigenvalue problem:
If ones multiplies each member of Equation (6) to the right by the transpose of the matrix of eigenvectors and considers normalized modes so that , the expression for the correlation matrix is given by:
The correlation matrix is then written as a sum of independent matrices defined as spatial autocorrelation patterns with proper modes as components. Since the eigenvectors form an orthonormal basis of the source space, can be written as:
where are the temporal modal amplitudes or the projection coefficients on the modal basis. is the i-th component of the k-th eigenvector ( ). Thus, coefficients are the root mean square of the acoustic pressure projected on the axis in the source space (i.e., the blade). According to Merces’ theory  , the projection coefficients form an orthogonal basis of the temporal space, and its root mean square corresponds to the eigenvalue . The eigenvalues represent the square of the sound pressure projected on the axis in the source space. The total sound pressure from the loading noise radiated to the receiver at time t is evaluated from Equations 1 and 8. It is given by:
Considering the orthogonal eigenvectors, the quadratic pressure at receiver is:
In this proper orthogonal decomposition, each mode does not contribute equally to the total quadratic pressure. The latter could then only be evaluated by taking into account the modes that contribute the most. This is reported by the accumulated acoustic energy of the first q modes, , divided by the total acoustic energy of all modes:
Once the modes that contribute the most are identified, summations in Equations (3) and (10) are limited up to instead of . Thus, Equations (9) and (10) become respectively:
In the previous equations, the spatial eigenvectors give information on the acoustic radiation of all dipole sources distributed over the surface S. Thus, considering Equation (13), the acoustic power defined by Equation (3) becomes:
where is the elementary surface of the ℓ-th receiver. In this study, is a constant. is the number of receivers, is the k-th eigenvector of the ℓ-th receiver, and is the i-th component of the eigenvector associated with the k-th eigenvalue .
2.3. Singular Value Decomposition (SVD) Approach
SVD investigates sound generation (i.e., the eigenvalues and the eigenvectors resulting from the POD) independent of the receiver . A global matrix gathering all the correlation matrices ( ) of the receivers is built. The global matrix of dimension is defined by:
Here, the receivers are distributed over a sphere of radius R around the source. The radius is large enough that the far field assumption is verified. The position of the receivers over the sphere is done according to ISO 3745  . SVD consists of decomposing the matrix in the following form:
where is the diagonal matrix of singular eigenvalues, and and are the matrices of left and right eigenvectors respectively. accounts for the difference of acoustic radiation of the sources in the free field of all receivers, and provides the average information on acoustic radiation of the sources in the free field of all receivers. One can show that and can be obtained by applying the POD to matrices and respectively. By retaining the first q energetic modes only, the expression of sound power (Equation (3)) based on the SVD can finally be expressed as:
where is the k-th eigenvalue of the matrix , and respectively represent the right and left i-th component eigenvectors, respectively associated with the eigenvalue .
In Figure 1, a centrifugal fan blade based on the modified NACA12 profile is used to illustrate our approach. The chord and wingspan equals and respectively. The blade is placed in a height periodic channel whose periodic faces are separated by the width . The cross section area of the channel equals , where and are the length and the width respectively. This area remains constant from the inlet until the
Figure 1. Geometry and boundary conditions.
outlet (Figure 1). The numerical domain comprises three volumes separated by interfaces: the inlet, the blade, and the outlet volumes, whose lengths equal , and respectively. The fluid is pure air and its physical properties are estimated at 25˚C ( and μ = 1.831 × 10−5 Pa・s). An airflow of is imposed on the inlet. The periodicity boundary conditions are applied to the periodic face and the no slip boundary conditions are imposed on the walls (blade, top and bottom). The purpose here is to investigate the possibility of both the POD and SVD approaches for identifying zones of the blade that are responsible for noise generation. Only frequencies in the 50 Hz to 10 kHz range are considered.
The calculation is performed on a hybrid spatial discretization realized with the free software Salome  . The discretization is composed of a structured and inhomogeneous mesh in the input and output volumes. In the blade volume and far from the blade, the mesh is unstructured .The mesh near the blade is structured with eight layers of geometrical progression 1.1 and total thickness . Around the blade, the dimensionless variables (flow direction), (spanwise direction), and (direction normal to the walls) are used. These values are sufficient according to Sagaut criteria  to correctly predict the boundary layer behavior when using LES method. This is also in accordance with de Villiers investigation  that employed a mesh satisfying these criteria and obtained numerical data close to the experiments. Recently, in one of our accepted articles for publication  , the influence of discretization on POD and SVD approaches has been studied. The coarse, medium and fine meshes describe a same noisiest area of the surface of the source. Thus, the mesh is composed of 1,100,000 cells, including cells on the blade.
The large eddy simulation method is used to simulate the internal flow channel with the free software OpenFOAM  . The one-equation subgrid model for turbulent kinetic energy is used  . The linear system is solved with the iterative preconditioned conjugate gradient method (PCG) with preconditioning DIC (diagonal incomplete Cholesky) for the pressure field, and the stable biconjugate gradient method (BiCGStab) with preconditioning DILU (diagonal incomplete LU) for the velocity field and turbulent kinetic energy. The backward second order scheme is used for the temporal resolution with a time step of for a Courant-Friedrichs-Lewy (CFL) number lower than 0.2. The gradient terms are calculated with a Gauss linear scheme, and the divergence terms are calculated with that of a Gauss vanLeerV 0.5 scheme for velocity and a Gauss upwind scheme for turbulent kinetic energy. Three weeks were necessary to reach convergence using 24 processors on Compute Canada-Sherbrooke  .
The second invariant of the velocity gradient tensor, named Q-criterion, was introduced by Hunt et al.  to better visualize the coherent structures of the flow characterized by a positive Q-criterion value. One observes in our configuration that there are no coherent structures in the region upstream the blade (Figure 2), despite the 5% intensity level imposed at the inlet. Coherent structures appear from a region located at 30% of the blade chord and develop gradually until the trailing edge. In this location, a strong intensification of the production of coherent structures that are then transported by the flow is particularly noticeable.
Once the calculation is converged, the wall pressure fluctuations on the blade are saved at time interval . samples of the pressure fluctuations are then the input data of both POD and SVD methods. The sound pressure, correlation matrix and modes are calculated for receivers placed on a sphere of radius R around the source, according to ISO 3745. The distance R is such that the far field hypothesis is verified (i.e. where is a characteristic length of the source and is the wave number based on the smallest cut-off frequency that is ). For this investigation, receivers were placed on the sphere of radius encompassing the source (Figure 3(a)). To better see the directivity of the radiation of the blade, 50 receivers were placed along a circle having the same radius as the
Figure 2. Isocontour Q-criterion from 0 to 1000: (a) All volume, and (b) Zoom on the middle volume and contour lines at the trailing edge of the blade.
Figure 3. Receivers positions: (a) Receivers on the sphere and (b) Receivers position on the circle.
sphere. The circle is in the plane perpendicular to the chord of the blade (i.e. in the plane . The 50 receivers are identified by their polar coordinates . The angular origin is from the axis. The positions of 4 of the 50 receivers are illustrated in Figure 3(b).
3.1. POD Analysis and Interpretation
The POD principle is to search an orthonormal basis of the m elementary sources of the discretized radiating blade’s surface. From this orthogonal basis, the q eigenvectors that contribute the most to the acoustic radiation at a given receiver are identified. The mapping of these q eigenvectors on the blade’s surface allows identifying its most radiant zones. Consequently, the application of the POD method on our stationary blade in the channel aims to minimize the number of modes necessary to understand the most radiant zones of the blade due to its interaction with the turbulent flow.
Indeed, Figure 4 presents the accumulated acoustic energy of the q first modes (Equation (11)) at receivers , , , and . One observes that the first 10 eigenvectors capture 99.60% of the total acoustic energy, whatever the receivers. As a result, these 10 eigenvectors could then be employed for the reconstruction of the total loading noise. This leads to a reduction in the number of modes (from to ) to be analyzed in order to identify the noisiest zones.
For example, Figure 5 shows the temporal evolution of the sound pressure of the loading noise reconstructed with (Equation (1)) and (Equation (9)), and (Equation (12)) at receivers , , , and . Acoustic pressures constructed with the first 10 modes (Equation (12)) of each receiver are in accordance with the reference acoustic pressure (Equation (1)). The maximum relative error equals 1.7% at receiver . The relative
Figure 4. POD percentage accumulated acoustic energy (eigenvalues, Equation (11)) of the first q modes divided by the total energy of all modes for receivers , , , and .
Figure 5. Comparison between total loading noise ( ) and loading noise reconstruction using the first 10 eigenvalues at receivers: (a) , (b) , (c) , and (d) .
error remains less than 0.1% for the other receivers.
One also observes that the temporal evolution of the sound pressure depends on the receiver location (Figure 5), which characterizes the directivity of radiation. It results in the POD being dependent on the receiver. For example, the maximum amplitude of sound pressure is almost equal to 0 for the receiver, and is equal to Pa for the receiver . The amplitude of the sound pressure signal is very small for receiver located at 0˚ in comparison with the others. In order to better see this signal a subfigure has been added on Figure 5(a).
Consider the receivers in Figure 3(b). We calculated the mean quadratic sound pressure with all modes (Equations (2) or (10)) and with the first 10 dominant modes (Equation (13)). First, One observes that the acoustic radiation is dipolar (Figure 6). The lowest sound pressure is noticed for receivers located at and , while the highest radiations are for receivers at and (Figure 6). This observation is in agreement with the directivity defined by Farassat  . According to this formulation (Equation (1)) the scalar product between the normal surface of the blade and the direction with the receiver is considered in the loading noise evaluation. Thus this product becomes null for the two former positions (i.e. egals 0˚ and 270˚) while equals unity for the two latter (i.e. egals 0˚ and 270˚).
In addition, the two curves of the mean quadratic sound pressure are similar, since the relative error remains less than 3% (Figure 6). This confirms the analysis of the reconstruction of sound pressure with the first 10 most important modes (Figure 5). When considering each term of the mean quadratic sound pressure (Equation (13)). We noted that the contribution of the modes does not necessarily appear in the ascending order (Figures 4-7. For example, considering receiver , the contribution of mode 1 is greater than the contribution of mode 0, even if the latter accumulates 71% of the energy, while mode 1
Figure 6. Directivity estimation in dB between mean quadratic sound pressure ( ) and mean quadratic sound pressure using the first ten eigenvalues at receivers of the Figure 3(b).
Figure 7. Mean quadratic sound pressure (product of eigenvalues and sum of the components of the eigenvector, Equation (13) for each of the first 15 modes of receivers , , , and .
only accumulates 23% (Figure 4). This contradiction is due to the phase shift of the components of an even eigenvector in calculating the mean quadratic sound pressure (Equation (13)). As a consequence the modes must be classified in ascending order when calculating either the mean quadratic sound pressure or the acoustic power. Thus, the reclassification of the modes was carried out for each receiver in Table 1.
One considers the reconstructed sound power (Equation (14)) with and . As shown in Table 2, the relative error of the acoustic power of the first five modes (Equation (14)) with respect to all modes (Equation (3)) is low and even lower with the first 10 modes (Equation (14)). This confirms that all useful acoustic information needed here is contained in the first 10 modes.
It is possible to project the eigenvectors that contribute the most to the acoustic power onto the blade surface, since the quadratic mean pressure is proportional to the sound power (Equations (13) and (14)). For example, when considering receivers , , , and , one observes that the POD method highlights the trailing edge as the region characterized by a great amplitude of the component , whatever the side of the blade (Figure 8 and Figure 9). The other regions of the blade do not seem to be affected by such amplitudes. Hence, mode 2 for receivers , 56˚, and 91˚, and mode 3 for receiver highly contribute to the acoustic power. The same observation is made for modes 3, 2, 4 and 9 for receivers , , , and respectively (last column of Figure 8 and Figure 9).
These previous observations are in agreement with the Q-criterion distribution (Figure 2), since according to this parameter, the turbulence is high near the trailing edge. The high velocity and pressure fluctuations in this region
Table 1. Classification of modes according to their contribution to the calculation of mean square acoustic pressure.
Table 2. Results of the sound power reconstructed with POD and relative error.
Figure 8. Mapping POD eigenvector on the overpressure side of the blade according to the classification given in Table 1: line 1, receiver ; line 2, receiver ; line 3, receiver ; line 4, receiver , and columns rank from most to least important modes from left to right.
would therefore make this region responsible for the loading noise.
It was shown here that both the radiated sound pressure and the sound power can be reconstructed with only a few POD modes. Instead of using all the modes
Figure 9. Mapping POD eigenvector on the underpressure side of the blade according to the classification given in Table 1: line 1, receiver ; line 2, receiver ; line 3, receiver ; line 4, receiver , and columns rank from most to least important modes from left to right.
proposed by the POD method (i.e. ), the first 10 dominant POD modes only were sufficient. In addition, the dipole character of the blade of radiation was noticed. The mapping of the major eigenvectors for several receivers demonstrated that the trailing edge is the region that radiates the most due to the high level of turbulence.
3.2. SVD Analysis and Interpretation
Like the POD, SVD methods based on CFD calculations and acoustic analogies were applied to identify the zones on a stationary blade in a channel that contribute the most to the sound power radiated in a subsonic regime. Using the SVD method, radiated sound power and pressure can be recovered using only the first few modes.
Figure 10 represents the acoustic energy that the q first modes accumulate divided by the total energy (Equation (11)). It emerges that 99.90% of the total acoustic energy is contained within the first 10 modes, which causes a reduction in the number of modes to be analyzed in order to understand the noisiest zones of the blade. The reconstructed sound power (Equation (17)) using the first 10 modes only is given in Table 3 and is compared with the sound power calculated
Figure 10. SVD percentage accumulated acoustic energy of the first modes divided by the total energy of all modes.
Table 3. Results of the sound power reconstructed with SVD and relative error.
using all the modes. The relative error equals 0.43%.
As shown previously for the POD method, the most energetic eigenvalues do not necessarily contribute the most to the reconstruction of the acoustic power. Modes 6, 10, 4, and 9 in order of their importance contribute the most (Figure 11).
In addition, the acoustic power estimated with only mode 6 equals (i.e., 50.6 dB) representing 93% of the total acoustic power. When adding the contribution of mode 10 to the contribution of mode 6, 96% of the total acoustic power is recovered. Thus, the number of modes to be analyzed to comprehend the most radiant zones of the blade regardless of receiver changes from modes to two modes (modes 6 and 10).
As for the POD method, it is possible to use singular value decomposition to map the right eigenvectors of the four main modes (i.e. modes 6, 10, 4, and 9) on the blade surface as shown in Figure 12. One notices that the amplitude for mode 6 is quite homogenous along the surface. As mentioned previously, this must be due to the turbulence that develops quite homogenously in the blade region; the entire surface therefore contributes equally to the loading noise. When considering the distribution for mode 10, one observes that the trailing edge is characterized by greater amplitudes, which is in accordance with the strong intensification of the turbulence in this region as shown by the
Figure 11. SVD acoustic power radiated by each mode.
Figure 12. Mapping SVD eigenvectors on the overpressure and underpressure side of the blade: line 1, overpressure side; line 2, underpressure side, and columns rank from most to least important modes.
Q-criterion (Figure 2). As a result, this region contributes largely to the radiated acoustic power through mode 10.
As a result, based on this observation, one could expect that the radiated sound power could be reduced by a proper modification of the blade geometry, or surface treatment, in these regions. This conclusion is yet to be validated in future work.
A CFD calculation was performed using OpenFOAM in order to estimate the wall pressure fluctuations on the surface of a blade located in a channel. The loading noise was then evaluated using the FW&H acoustic analogy and decomposed using both the POD and SVD methods. It was observed that the sound power and the radiated pressure were recovered using these two methods, even if only a few modes were considered. In the configuration studied here, 10 POD modes or 1 SVD mode can be sufficient to predict the radiated sound pressure. Whatever the approach, the trailing edge of the blade is distinguished. An acoustic treatment or a geometrical modification of this zone can therefore influence the blade’s acoustic radiation.
It was shown that the POD approach provides information about the directivity of the source, which is not the case when using the SVD approach. However, the latter approach would be employed when a global noise reduction is sought, no matter the direction of the sound. On the contrary, if a decrease is preferred in a particular direction, the POD method should be used.
This work was supported by Natural Sciences and Engineering Research Council of Canada (NSERC). The authors also wish to thank Compute Canada in Sherbrooke, Quebec, for its assistance.