*2.2. Validation of LES for Low Mach Number Turbulent Flows around NACA0012 Airfoil*

The calculation target of this study is a three-dimensional flow around a two-dimensional wing of the NACA0012 airfoil, which is a typical streamlined object. NACA0012 airfoil is defined by the following equation and its geometry is given in Figure 1.

$$\pm y/\mathbb{C} = 0.6 \times \left[ 0.2969 \sqrt{\mathbf{x}/\mathbb{C}} - 0.1260 (\mathbf{x}/\mathbb{C}) - 0.3516 (\mathbf{x}/\mathbb{C})^2 + 0.2843 (\mathbf{x}/\mathbb{C})^3 - 0.01015 (\mathbf{x}/\mathbb{C})^4 \right]. \tag{18}$$

**Figure 1.** The profile of NACA0012 airfoil.

In order to validate our method, we conducted a LES of the compressible flow around NACA0012 airfoil with the computational setup matching the experimental setup of Miyazawa et al. [31], in which the Mach number was 8.75 × <sup>10</sup><sup>−</sup>2, the angel of attack was 9◦, and the Reynolds number was 2 × <sup>10</sup>5, which was based on the uniform velocity in the streamwise direction and the chord length. Figure 2 shows the computational domain and computational boundary conditions in the present study. For the Cartesian coordinate system, *x*, *z*, and *y* were defined in the streamwise direction, spanwise direction, and vertical direction, respectively. We used a C-type boundary-fitted grid in the *x*–*y* plane. Thus, a general curvilinear coordinate system had to be applied and represented as *ξ*, *η*, *ζ* for all computations, in which the direction following the mainstream surface of the airfoil is denoted as *ξ*, the direction away from the surface of airfoil is defined as *η*, and the spanwise direction is represented as *ζ*. The computational domain size is defined as the following: the length of the wake, the semidiameter of the C-type grid, and the spanwise extent were 11*C*, 11*C*, and 0.5%*C*, respectively. Here, *C* is chord length. For the mesh, there were 1600 grid points in the *ξ* direction, in which 800 was on the airfoil surface, 160 grid points were in the *η* direction, and 60 grid points were in the *ζ* direction. As shown in Figure 2, the non-slip boundary condition was applied on the airfoil surface. The inflow was the uniform stream without disturbance. The convective outflow was applied for the outflow boundary condition. The spanwise direction boundary condition is defined as the periodic. For the variables in the *η* direction, zero gradient was employed at the op and bottom boundary. Particular attention was paid for the non-reflective boundary condition by [32] at the far boundary.

A second-order central finite-difference discretization scheme was applied in the equation of motion for the diffusion terms. In order to improve the numerical stability in high Reynolds number flow, a QUICK method was applied for the convective terms in the general curvilinear coordinate system. For coupling the pressure field and the continuity equation, the fractional method was applied. The second-order Adams–Bashforth method is used to the convective term and the viscous term for the advancement method of time in the equation of motion. For the dynamic procedure, the test filter was used in the streamwise direction and spanwise direction with second-order accuracy and the test-to-grid filter ratio Δ˜ /Δˆ = 2. The present numerical method and computer program have been tested extensively in several turbulent flows [33–37].

**Figure 2.** Computational domain and boundary conditions.

All simulations were computed on an NEC SX-8R supercomputer of Cybermedia Center, Osaka University with the time step *dt* = 0.0495*H*/*Uc*. The greater part of the total effort toward calculations was spent on solving the pressure equation through the residual cutting method [38]. All simulations were run until the flow fields that were fully developed and the first-order, and second-order statistics exhibited adequate convergence. The time-averaging, and, the spatial averaging in the spanwise direction were used for collecting the data. Figures 3 and 4 show the distribution of the mean pressure coefficient *Cp* on the airfoil surface and the average pressure coefficient fluctuation *Cp rms* on the suction side of the airfoil, in which *Cp* is defined by the freestream pressure *p*0.

$$\mathcal{C}\_p = \frac{p - p\_0}{\frac{1}{2}\rho l I\_0^2}.\tag{19}$$

Then, the root-mean-square of *Cp* is defined as *Cp rms*. In order to demonstrate the validity of our model, the results of our LES are compared with the measurement by [31]. In Figures 3 and 4, our computational results are in good agreement with the experimental data.

**Figure 3.** Mean pressure coefficient *Cp*.

**Figure 4.** Mean fluctuation of *Cp*.

#### **3. Derivation of a New Acoustic Model**

A new acoustic model was derived by using the compressible continuity and Navier–Stokes equations of the flow field; i.e.,

$$\frac{\partial \rho}{\partial t} + \frac{1}{J} \frac{\partial}{\partial \xi^k} \left( J \frac{\partial \xi^k}{\partial x\_i} \rho u\_j \right) = 0,\tag{20}$$

$$\frac{\partial \mu\_i}{\partial t} + \frac{1}{J} \frac{\partial}{\partial \xi^k} \left( J \frac{\partial \xi^k}{\partial x\_j} \rho u\_j u\_i \right) - \frac{\partial \xi^k}{\partial x\_i} \frac{\partial p}{\partial \xi^k} = \frac{\partial \xi^k}{\partial x\_j} \frac{\partial \sigma\_{i\bar{j}}}{\partial \xi^{k\bar{j}}} \tag{21}$$

$$\sigma\_{ij} = 2\mu \left( S\_{ij} - \frac{1}{3} \delta\_{ij} S\_{kk} \right), \quad S\_{ij} = \frac{1}{2} \left( \frac{\partial \xi^k}{\partial \mathbf{x}\_i} \frac{\partial u\_j}{\partial \xi^k} + \frac{\partial \xi^k}{\partial \mathbf{x}\_j} \frac{\partial u\_i}{\partial \xi^k} \right). \tag{22}$$

Variables *ρ* and *p* are decomposed as

$$
\rho = \rho\_0 + p', \quad p = p\_0 + p', \tag{23}
$$

where (0) indicates the mean component of variables; (- ) means the perturbation component of variables, due to in low Mach number turbulent flow

$$p\_0 \gg p \,, \quad p\_0 \gg p \,. \tag{24}$$

Rewriting Equations (20) and (21) by using Equation (23) leads to

$$\frac{\partial \rho^{\prime}}{\partial t} + \frac{1}{J} \frac{\partial}{\partial \tilde{\xi}^{k}} \left( J \frac{\partial \tilde{\xi}^{k}}{\partial \mathbf{x}\_{i}} \rho^{\prime} \mathbf{u}\_{j} \right) = -\rho\_{0} \frac{\partial \tilde{\xi}^{k}}{\partial \mathbf{x}\_{j}} \frac{\partial \mathbf{u}\_{j}}{\partial \tilde{\xi}^{k}} \tag{25}$$

$$\frac{\partial u\_i}{\partial t} + \frac{1}{J} \frac{\partial}{\partial \xi^k} \left( J \frac{\partial \xi^k}{\partial x\_j} u\_j u\_i \right) - \frac{1}{\rho} \frac{\partial \xi^k}{\partial x\_i} \frac{\partial p^r}{\partial \xi^k} = \frac{1}{\rho} \frac{\partial \xi^k}{\partial x\_j} \frac{\partial \sigma\_{i\bar{j}}}{\partial \xi^k}. \tag{26}$$

After multiplying *ui* to Equation (25) and *ρ*, to Equation (26), Equations (25) and (26) are added. Namely,

$$\frac{\partial \rho \, u\_i}{\partial t} + \frac{1}{J} \frac{\partial}{\partial \xi^k} \left( J \frac{\partial \xi^k}{\partial x\_i} \rho \, u\_i u\_j \right) + \frac{\rho}{\rho} \frac{\partial \xi^k}{\partial x\_i} \frac{\partial \rho^r}{\partial \xi^k} = \frac{\rho}{\rho} \frac{\partial \xi^k}{\partial x\_j} \frac{\partial \sigma\_{ij}}{\partial \xi^k} - \rho\_0 u\_i \frac{\partial \xi^k}{\partial x\_j} \frac{\partial u\_j}{\partial \xi^k}. \tag{27}$$

Taking the divergence of Equation (27), subtracting *<sup>∂</sup> <sup>∂</sup><sup>t</sup>* of Equation (25), and then adding −*c*<sup>2</sup> 0 1 *J ∂ ∂ξk* - *<sup>γ</sup>kl ∂ρ*, *∂ξl* concerning two sides of the equation obtained, the wave equation with source terms is finally obtained as the following

*∂*2*ρ*, *<sup>∂</sup>t*<sup>2</sup> <sup>−</sup> *<sup>c</sup>*<sup>2</sup> 0 1 *J ∂ ∂ξ<sup>k</sup> <sup>γ</sup>kl ∂ρ*, *∂ξ<sup>l</sup>* <sup>=</sup> <sup>1</sup> *J ∂ ∂ξ<sup>k</sup> <sup>γ</sup>kl <sup>∂</sup> ∂ξ<sup>l</sup> ρ*, *uiuj* <sup>−</sup> *<sup>c</sup>*<sup>2</sup> 0 1 *J ∂ ∂ξ<sup>k</sup> <sup>γ</sup>kl ∂ρ*, *∂ξ<sup>l</sup>* <sup>+</sup> *<sup>ρ</sup>*, *ρ* 1 *J ∂ ∂ξ<sup>k</sup> <sup>γ</sup>kl ∂ρ*, *∂ξ<sup>l</sup>* <sup>+</sup> *∂ξ<sup>l</sup> ∂xi ∂p*, *∂ξ<sup>l</sup> ∂ξ<sup>k</sup> ∂xi ∂ ∂ξ<sup>k</sup> ρ*, *<sup>ρ</sup>* <sup>−</sup> *<sup>ρ</sup>*<sup>0</sup> *∂ ∂t ξk ∂xi ∂ui ∂ξ<sup>k</sup>* − *ρ*, *ρ* 1 *J ∂ ∂ξ<sup>k</sup> <sup>γ</sup>kl <sup>∂</sup> ∂ξ<sup>l</sup> σij* (28) <sup>−</sup> *∂ξ<sup>k</sup> ∂xi ∂ξ<sup>l</sup> ∂xj ∂σij ∂ξ<sup>k</sup> ∂ ∂ξ<sup>l</sup> ρ*, *ρ* + *ρ*<sup>0</sup> *∂ξ<sup>l</sup> ∂xi ∂ξ<sup>k</sup> ∂xj ∂ui ∂ξ<sup>l</sup> ∂uj ∂ξ<sup>k</sup>* <sup>+</sup> *<sup>ρ</sup>*0*ui* 1 *J ∂ ∂ξ<sup>k</sup> <sup>γ</sup>kl <sup>∂</sup>uj ∂ξ<sup>l</sup>* .

The first term and the last term of the right-hand side of Equation (29) can be rewritten as

$$\frac{1}{\mathcal{I}}\frac{\partial}{\partial\widetilde{\xi}^{\mathbb{K}}}\left(\gamma^{kl}\frac{\partial}{\partial\widetilde{\xi}^{l}}\rho\cdot\mu\_{i}\mu\_{j}\right) + \rho\_{0}\mu\_{l}\frac{1}{\mathcal{I}}\frac{\partial}{\partial\widetilde{\xi}^{\mathbb{K}}}\left(\gamma^{kl}\frac{\partial\mu\_{l}}{\partial\widetilde{\xi}^{l}}\right) = \frac{1}{\mathcal{I}}\frac{\partial}{\partial\widetilde{\xi}^{\mathbb{K}}}\left(\gamma^{kl}\frac{\partial}{\partial\widetilde{\xi}^{l}}\rho\mu\_{i}\mu\_{j}\right) + \rho\_{0}\mu\_{l}\frac{\partial\widetilde{\xi}^{\mathbb{K}}}{\partial\widetilde{x}\_{l}}\frac{\partial}{\partial\widetilde{\xi}^{l}}\left(\frac{\partial\widetilde{\xi}^{l}}{\partial\widetilde{x}\_{l}}\frac{\partial\mu\_{l}}{\partial\widetilde{\xi}^{l}}\right). \tag{29}$$

Upon substituting Equation (29) into Equation (29), the wave equation with a lump of source terms is obtained as follows

$$
\begin{split}
\frac{\partial^2 \rho^\prime}{\partial t^2} - \varepsilon\_0^2 \frac{1}{\overline{l}} \frac{\partial}{\partial \xi^k} \left(\gamma^{kl} \frac{\partial \rho^\prime}{\partial \xi^l}\right) &=& \frac{1}{\overline{l}} \frac{\partial}{\partial \xi^k} \left(\gamma^{kl} \frac{\partial}{\partial \xi^l} \rho u\_i u\_j\right) - \varepsilon\_0^2 \frac{1}{\overline{l}} \frac{\partial}{\partial \xi^k} \left(\gamma^{kl} \frac{\partial \rho^\prime}{\partial \xi^l}\right) + \frac{\rho^\prime}{\rho} \frac{1}{\overline{l}} \frac{\partial}{\partial \xi^k} \left(\gamma^{kl} \frac{\partial \rho^\prime}{\partial \xi^l}\right) \\
&+& \frac{\partial^2\_\xi}{\partial x\_i} \frac{\partial \rho^\prime}{\partial \xi^l} \frac{\partial \xi^k}{\partial x\_i} \frac{\partial}{\partial \xi^k} \frac{\rho^\prime}{\rho} - \frac{\rho^\prime}{\rho} \frac{1}{l} \frac{\partial}{\partial \xi^k} \left(\gamma^{kl} \frac{\partial}{\partial \xi^l} c\_{ij}\right) \\
&- \frac{\partial^2\_\xi}{\partial x\_i} \frac{\partial \xi^l}{\partial x\_j} \frac{\partial \gamma^l}{\partial \xi^k} \frac{\partial \gamma^l}{\partial \xi^l} \left(\frac{\rho^\prime}{\rho}\right) - \rho\_0 \frac{D}{Dt} \left(\nabla \cdot \mathbf{u}\right) + \rho\_0 \left(\nabla \cdot \mathbf{u}\right)^2,\end{split}
\tag{30}
$$

where *<sup>D</sup> Dt* <sup>=</sup> *<sup>∂</sup> <sup>∂</sup><sup>t</sup>* + *uj ∂ <sup>∂</sup>xj* means the material derivative. Here, assuming a low Mach number flow (*ρ*<sup>0</sup> *<sup>ρ</sup>*, ) and a high Reynolds number flow (*δij* 1 ), a wave equation with source terms can be finally obtained as

$$\frac{\partial^2 \rho}{\partial t^2} - c\_0^2 \frac{1}{\hbar^3} \frac{\partial}{\partial \dot{\xi}^l} \left( \gamma^{kl} \frac{\partial \rho}{\partial \dot{\xi}^l} \right) = \frac{1}{\hbar} \frac{\partial}{\partial \dot{\xi}^k} \left( \gamma^{kl} \frac{\partial}{\partial \dot{\xi}^l} \rho u\_l u\_j \right) - c\_0^2 \frac{1}{\hbar} \frac{\partial}{\partial \dot{\xi}^k} \left( \gamma^{kl} \frac{\partial \nu}{\partial \dot{\xi}^l} \right) - \rho\_0 \frac{D}{Dt} \left( \nabla \cdot \mathbf{u} \right) + \rho\_0 \left( \nabla \cdot \mathbf{u} \right)^2. \tag{31}$$

The form of our wave Equation (31) is similar to the form of the following Lighthill's equation [11]. Lighthill's equation is expressed as

$$\frac{\partial^2 \rho^\prime}{\partial t^2} - c\_0^2 \frac{1}{I} \frac{\partial}{\partial \xi^{kl}} \left( \gamma^{kl} \frac{\partial \rho^\prime}{\partial \xi^{l\prime}} \right) = \frac{1}{I} \frac{\partial}{\partial \xi^{kl}} \left( \gamma^{kl} \frac{\partial T\_{ij}}{\partial \xi^{l\prime}} \right), \tag{32}$$

where

$$T\_{ij} = \rho u\_i u\_j + \delta\_{ij} \left[ (p - p\_0) - c\_0^2 (\rho - \rho\_0) \right] - \sigma\_{ij}.\tag{33}$$

Equation (31) and Lighthill's Equation (32) both seem to be wave equations having the source terms on the right hand side. But there are unknowns in *Tij* for Lighthill's Equation (32). Especially in case of low Mach number flows, the simplification of *Tij* ≈ *ρ*0*uiuj* works well. Hence the indirect method, which solves the acoustic field using *Tij* given by the incompressible solver, becomes reasonable. On the other hand, from our wave Equation (31), the utilization of ∇ · **u** for the sound source is possible. Actually, the divergence of the velocity vector plays an important role in the flow field even in the low Mach number [34]. Therefore, we rearranged the acoustic wave equation by exposing ∇ · **u** intentionally and then obtained Equation (31), so that the indirect method could be extended from the zero Mach number region to low or moderate Mach number region.

For easily applying to engineering practices in the industry, such as large-scale wind turbine blades and underlining the key role of ∇ · **u** in the flow field, our acoustic model Equation (31) can be further approximated as

$$\frac{\partial^2 \rho^\prime}{\partial t^2} - c\_0^2 \frac{1}{J} \frac{\partial}{\partial \xi^k} \left( \gamma^{kl} \frac{\partial \rho^\prime}{\partial \xi^l} \right) \approx -\rho\_0 \frac{D}{Dt} \left( \nabla \cdot \mathbf{u} \right) + \rho\_0 \left( \nabla \cdot \mathbf{u} \right)^2. \tag{34}$$

The Lighthill's Equation (32) is exact, but some kind of turbulence modeling for computing the sound source is required when the Reynolds number is high. Especially in case of low Mach number flows, the simplification of *Tij* ≈ *ρ*0*uiuj* is utilized so that the solver becomes incompressible. Powell's sound source model [19] is derived under the assumption of incompressible flows. On the other hand, although some approximation is used in the derivation of the acoustic equation, our sound source model (34) does not need any modeling process. Furthermore, it can consider the influence of the compressibility effect. Thus, the behavior of sound source due to the variation of Mach number is expected to be reproduced appropriately.

Applying our acoustic theory and the compact green function [22], the sound pressure in the far-field is solved from the wave Equation (34). That is, assuming that the observation point *x* is sufficiently far from the sound source area of the objet *y*, the object keep stationary and the velocity of the surface *S* of the object is zero, the wave Equation (34) is able to be solved applying the compact green function:

$$p\_a(\mathbf{x}, t) = \rho\_0 \iint \frac{D}{Dt} (\nabla \cdot \mathbf{u}) \cdot \nabla G(\mathbf{x}, y, t - \mathbf{r}) d^3 y d\mathbf{r},\tag{35}$$

where *pa* denotes the pressure sound, *G* means the compact green function and is represented as

$$G(\mathbf{x}, y, t - \tau) = \frac{1}{4\pi|\mathbf{x}|} \delta\left(t - \tau \frac{|\mathbf{x}|}{c\_0}\right) + \frac{\mathbf{x} \cdot \mathbf{y}}{4\pi c\_0 |\mathbf{x}|^2} \frac{\partial}{\partial t} \delta\left(t - \tau - \frac{|\mathbf{x}|}{c\_0}\right). \tag{36}$$

Here, *<sup>Y</sup>* means Kirchhoff vector and is calculated through ∇2*Y<sup>i</sup>* = 0. Then substitution of Equation (36) into Equation (35) results in

$$p\_{\mathbf{a}}(\mathbf{x},t) = \frac{\rho\_0 \mathbf{x}\_i}{4\pi \varepsilon c\_0 |\mathbf{x}|^2} \frac{\partial}{\partial t} \int \frac{D}{Dt} (\nabla \cdot \mathbf{u}) \left(\mathbf{y}, t - \frac{|\mathbf{x}|}{c\_0}\right) \nabla Y\_i(\mathbf{y}) d^3 y. \tag{37}$$

From the Equation (37), the sound pressure *p*, *<sup>a</sup>* can be obtained from conducting the Fourier transform. Finally the sound pressure level (SPL) can be calculated from the following equation through applying the solved *p*, *a*:

$$SPL(dB) = 10 \times \log \left(\frac{p\_a^2}{P\_b^2}\right) \,\text{}\,\tag{38}$$

where *Pb* means the reference sound pressure. In general, the magnitude of *Pb* is *Pb* = <sup>2</sup> × <sup>10</sup>−5*Pa*.

#### **4. Results of the Acoustic Field around NACA0012 Airfoil**
