Research on flow-sound separation algorithm of aerodynamic noise based on immersed boundary method

2.1. Control equations

In this paper, the hybrid method which was based on hydrodynamic force and sound separation [25, 26] was adopted to directly calculate aerodynamic noise under low Mach number. All flow variables were decomposed into the incompressible parts and compressible perturbation parts in this method:

Control equations of the incompressible part was N-S equation (INS), which represented flow field information and sound pulsation. And the other compressible effects were expressed by perturbation variable ($\text{'}$). N-S equation was written as follow:

Sound perturbation variables can be acquired by solving the linearized compressible perturbation equation (LPCE) after obtaining the incompressible flow solution [27]. The vector of LPCE equation can be written as follow:

INS/LPCE hybrid method has been successfully verified by the basic acoustic problems such as dipole and quadrupole [27]. Meanwhile, it was also applicable to turbulent flow noise problem [7, 9]. The left side of LPCE equation indicated the propagation and emission of sound waves, while its right side only included sound source item.

2.2. Numerical method

The incompressible N-S Eq. (3) was solved by several steps. The spatial derivatives were handled by two order central difference [28]. Pressure Poisson equation was solved by means of multiple grid method based on Line-Gauss-Seidel (LGS) matrix. LPCE spatial discretization was subject to six-order compact central difference [31], and time integral was subject to four-order Runge-Kutta. Three-order and four-order boundary were adopted near the immersed solid boundary [31]. The implicit spatial filtering [32] was adopted to restrain the high frequent errors, because there was no dissipation error in compact central difference. Ten-order filtering was adopted in the internal areas in the paper, and the filtering near the boundaries reduced to two-order from eight-order successively. The compact central difference and filtering format of implicit spatial can be solved by means of a tridiagonal matrix.

2.3. Immersed boundary method

The flow problems of the complicated immersed boundaries were solved by the incompressible N-S equation [28]. In this method, the surface of the immersed object was indicated by the unstructured surface grid of triangular elements. Before integral for control equations, all the elements inside the object were marked as “body points”, and the other points outside the object were marked as “fluid points”. Any element which was near fluid elements was marked as “image point”, as shown in Fig. 1(a). Wall boundary conditions were realized by assigning proper values for the image points, as shown in Fig. 1(b).

Fig. 1a) Diagram of image points and b) diagram of boundary conditions

a)

b)

The higher order immersed boundary method [23] adopted in sound solution was high order polynomial interpolation which was based on weighted Least Squares Error. Image point values determined whether it met the boundary conditions at Body Intercept (BI) points. In particular, a general variable $\varphi$ which was near BI points $({\mathbf{x}}_{BI},{\mathbf{y}}_{BI},{\mathbf{z}}_{BI})$was estimated by using $N$-order polynomial:

where, $x\mathrm{\text{'}}=x-x_{BI}$, $y\mathrm{\text{'}}=y-y_{BI}$, $z\mathrm{\text{'}}=z-z_{BI}$, $c_{ijk}$ was an unknown coefficient and was expressed as follows:

Sizes of coefficients for three-order polynomial ($N=$ 3) was 10 and 20 under two-dimensional and three-dimensional conditions, respectively. Coefficients of polynomials with different orders were described in reference [23]. To determine these coefficients, the data of the flow field points near BI points needed to be acquired. According to method of Luo [30], the flow field values in the circle or ball (under three-dimensional conditions) with the object point as the center and the radius of $R$ can be selected as the data. For $M$ data points, coefficient $c_{ijk}$ was determined by the minimum weighted error estimation:

where, $({\mathbf{x}}_m^{\text{'}},{\mathbf{y}}_m^{\text{'}},{\mathbf{z}}_m^{\text{'}})$ was the $m$th data point, $w_m^{}$ was weight function. According to experience of the previous researchers [30], cosine weight function was used in this paper. To make the least square problem appropriate, radius $R$ should be selected in an appropriate way which can ensure that there were more sizes of data points than that of coefficients. To find the values of image points and BI points, the first data point was a image point, the other ($M-$1) data points shall be found in the fluid region (as shown in Fig. 1(b)). The precise solution of least square problem (Eq. (9)) was as follow:

where the superscript “+” indicated the generalized inverse of matrix, vector $\mathbf{c}$ and $\varphi$ included coefficient $c_{ijk}$ and data $\varphi \left({\mathbf{x}}_m^{\text{'}},{\mathbf{y}}_m^{\text{'}},{\mathbf{z}}_m^{\text{'}}\right)$, respectively. $\mathbf{W}$ and $\mathbf{V}$ indicated weight matrix and Vandermonde matrix, respectively. $\left({\mathbf{x}}_1^{\text{'}},{\mathbf{y}}_1^{\text{'}},{\mathbf{z}}_1^{\text{'}}\right)$ indicated the first image point. Coefficient $c_{ijk}$ can be written as the linear combination of $\varphi \left({\mathbf{x}}_m^{\text{'}},{\mathbf{y}}_m^{\text{'}},{\mathbf{z}}_m^{\text{'}}\right)$ after solving Eq. (10). According to Eq. (8), coefficient $c_{ijk}$ indicated the value and derivative value of $\left({\mathbf{x}}_{BI},{\mathbf{y}}_{BI},{\mathbf{z}}_{BI}\right)$:

Therefore, for Dirichlet and Neumann boundary conditions given on wall surface, Image point values can be acquired by Eq. (10) and Eq. (11). The details of the equation for immersed boundary conditions can refer to reference [23].

2.4. Disposing the freshly cleared point

Disposing the moving object on the fixed grid would generate “freshly cleared points”, as shown in Fig. 2(a). As the time-dependent variable values of the freshly cleared points did not need to participate in integral of control equations, they can be obtained through interpolation among the near elements [33]. In solving the incompressible flow, the variable values of the new time layer were acquired by iteration of the momentum equation [28]. In solving acoustic solution, the values of the freshly cleared points were acquired through the high order multinomial Eq. (7). The entire assignment process was similar to the disposal of image points in Section 2.3. However, for freshly cleared points, the center of the data points was the freshly cleared point $(x_{FC},y_{FC},z_{FC})$, where $x\mathrm{\text{'}}=x-x_{FC}$, $y\mathrm{\text{'}}=y-y_{FC}$, $z\mathrm{\text{'}}=z-z_{FC}$, as shown in Fig. 2(b). To avoid iteration, only the non-updated elements and fluid elements are selected as the data points to minimize Least Squares Error. Solving Eq. (10) will acquire the coefficient of the approaching multinomial, and then the value of the updated element can be directly given by the first coefficient, for instance:

3.1. Turbulent flow noise of two tandem circular cylinders

As shown in Fig. 3, the mentioned method was applied to calculating the flow-induced noise of two tandem circular cylinders. This problem was a standard calculation example for researching airframe noise [34]. The inflow velocity $U_0=\mathrm{}$44 m/s, and the corresponding Reynolds number $\mathbf{R}{\mathbf{e}}_D=$1.66×10⁵. Mach number $M=$0.128 was suitable for the hybrid algorithm. During computing noise by means of this method, the computational domain size was 30$D\le \mathbf{x}\le$40$D$, –40$D\le \mathbf{y}\le$40$D$. And span-wise length ${\mathbf{L}}_z=\text{3}D$. The periodic boundary condition was applied in the span-wise direction. Sizes of the non-homogeneous Cartesian grids were 768×384×32. Fig. 4 showed the computational result in $\mathbf{x}$-$\mathbf{y}$ plane. The periodic boundary conditions were not marked in the figure, and it can be easily understood. We cared more about the cylinder wake vortex shedding and its mutual influences, and the periodic boundary conditions were applied to up, down, left and right planes in the computational domain. This was a very conventional method used in three-dimensional computation:

where $\mathrm{\Omega }$ and $S$ were the vortex and strain rate tensor, respectively. Under the current Reynolds number, the distance between the two circular cylinders was $\mathbf{S}=$ 3.7$D$. The wake flow of the upstream circular cylinder was in rolling up status before arriving at the downstream circular cylinder. The shedding vortex of the upstream circular cylinder and that of the downstream circular cylinder interacted with each other. The characteristics of the entire flow field were similar to that of the higher Reynolds number [35]. The change curve of aerodynamic coefficient with time was shown in Fig. 5. The coefficient average value and mean square root (rms) was shown in Table 1. It was shown that the aerodynamic force fluctuation of the downstream circular cylinder was more serious, because of the interference of the shedding vortex of the upstream circular cylinder. The main frequency of vortex shedding was $\mathit{S}\mathit{t}=$0.196. The aerodynamic force was higher than that of the experimental result under the higher Reynolds number ($\mathbf{R}{\mathbf{e}}_D=$1.66×10⁵) [34-36]. Although the shedding frequency ($\mathit{S}\mathit{t}=$0.196) was lower than the experimental result of NASA [36-38] ($\mathit{S}\mathit{t}=$0.234), it was very close to the direct numerical simulation result ($\mathit{S}\mathit{t}=\mathit{}$0.18, $\mathbf{R}{\mathbf{e}}_D=$ 1000) of Papaioannou [37] and the experimental result ($\mathit{S}\mathit{t}=$ 0.19, $\mathbf{R}{\mathbf{e}}_D=$ 4000) of Igarashi.

Fig. 3Geometrical model of two tandem circular cylinders

Fig. 4Vortex structure of two tandem circular cylinder

Fig. 5Change curve of aerodynamic coefficient with time

a) Upstream circular cylinder

b) Downstream circular cylinder

Computation of sound field was based on the incompressible flow, with LPCE as control equations. The span-wise was assumed to be periodic during computing the flow field. As scale of the sound wave length was larger than flow scale under current Mach number ($M=\mathrm{}$0.128), span-wise length ${\mathbf{L}}_z=\text{3}D$ was not enough to computer three-dimensional sound field. Based on the research of Seo [7], three-dimensional sound field was directly related to the central plane in span-wise direction, as a result, the two-dimensional sound field at zero span-wise direction (${\mathbf{k}}_z=\text{0}$) was computed. The computational result was modified to three-dimension by Oberai formula [39]. The contour line of time-average expansion rate for sound field was shown in Fig. 6(a). It can be observed that the wave-length of the main frequency sound wave was approximate 40.5$D$, and the high frequency components caused by turbulent fluctuation can also be observed. The sound pressures were observed at three points, which were $A$(–8.33$D$, 27.815$D$), $B$(9.11$D$, 32.49$D$) and $C$(26.55$D$, 27.815$D$), respectively. The positions of the observation points were the same as that of the sensors arranged in NASA experiment. The power spectral density of the sound pressure pulsation at each observation point was shown in Fig. 6(b). The spectrogram was promoted to the central plane data of three-dimensional result through modification [7]. Sound pressure frequency spectrum was constituted of broadband noise and pure tones with large amplitudes at main frequency. It was consistent with the experimental data [36]. The computational domain in $\mathbf{x}$-$\mathbf{y}$ plane of sound field was kept consistent with that of the flow field, the grid size was 500×400. The grid resolution near sound field was as twice as that of the flow field. It was more accurate at far field to accurately distinguish the high frequency sound waves.

Fig. 6a) Contour line of the instantaneous expansion rate of sound field and b) the power spectral density (PSD) of sound pressure at three observation points

a)

b)

Fig. 7The power spectral density of three observation points at actual span length (16D) as modified, the solid line was the computational result, and the dot dash line with symbols was NASA experimental data [34, 36]

This paper attempted to compare the sound results with the experimental one. Because sound prediction was carried out under small span-wise length (${\mathbf{L}}_z=\text{3}D$), so it must be modified to the actual span-wise length ($\mathbf{L}=\text{16}D$) to be compared with the experimental value. This required the information related to span-wise length. In this paper, the span-wise length data provided by the experiments [34] was adopted. It was found that the length related to span-wise was longer than the span-wise length only under main shedding frequency. Based on the modification formula of Seo, the modification value of main shedding frequency was $a+$ 9.4 [dB], while the other frequencies were $a+$ 7.2 [dB]. The modified power spectral density and experimental result were shown in Fig. 7. Because Reynolds number of the computation was different from the experimental one, the frequency spectrum was not identical well, in particularly the wave peak and total amplitude. However, it can be found that the change trend was consistent each other. For instance, there were second-order harmonic and third-order harmonic at point A simultaneously. But third-order harmonic was obvious at point B, while second-order harmonic was obvious at point C.

3.2. Flow noise of the rudimentary landing gear

The rudimentary landing gear (RLG) model used in the present study had a modular design as shown in Fig. 8(a). The model had three main components, the vertical post, the truck, and the wheels. The model had one instrumented indexed wheel, one instrumented non-indexed wheel, and two dummy wheels. The fabrication of the landing-gear model was carried out using composites. The post was fabricated using carbon-fiber-reinforced plastics, while the truck and the wheels were fabricated using glass-fiber-reinforced plastics. Using composites for fabrication ensured that the model weight was minimal without compromising the strength of the model. The total grid number of RLG model was 187742 triangular elements.

Fig. 8a) Geometric model of the rudimentary landing gear and b) Computational domain grid in x-y plane

a)

b)

Fig. 9 provided the topology and dimensions of the rudimentary landing gear model normalized with the wheel diameter $D$. Wheel width was 0.37$D$ and shoulder radius was 0.115$D$. Four-wheel truck wheelbase was 1.16$D$ and track was 0.88$D$. Transverse axles square with side was 0.3$D$. Longitudinal beam was 0.3$D$ high and 0.25$D$ wide. The vertical post was 0.25$D$ square.

To verify the applicability of this method for complicated profiles, the flow noise of the rudimentary landing gear was calculated in this section. The computational domain was a rectangle which sizes was $\text{0}D\le \mathbf{x}\le \text{12}D$, $\text{0}D\le \mathbf{y}\le \text{6}D$, $\text{0}D\le \mathbf{z}\le \text{5}D$ ($D$ was the diameter of the wheel). The total number of non-uniform Cartesian grid was 512×256×256. The computational domain grid in $\mathbf{x}$-$\mathbf{y}$ plane was shown in Fig. 8(b). $\mathbf{R}{\mathbf{e}}_D=$ 2000 and $M=\mathrm{}$0.3. Fig. 10(a) was the instantaneous vortex structure of the rudimentary landing gear. It was found that there was a complicated three-dimensional vortex structure in the wake flow area. Contour map of total pressure pulsations in different planes was shown in Fig. 10(b).

Based on flow field characteristics in Fig. 10, the sound field around the rudimentary landing gear was computed. Sound pressure levels at three observation points in Fig. 11 were extracted, as shown in Fig. 12.

Fig. 9Model topology and dimensions

Fig. 10Instantaneous flow field and sound field

a) Vortex structure of the flow field

b) Contour map of total pressure pulsation

Fig. 11Positions of three observation points

Fig. 12Sound pressure levels of three observation points

According to Fig. 12, sound pressure levels of observation point 2 was obviously higher than that of other two points. Its reasons were mainly as shown in the following. According to Fig. 10(a), the flow field of observation point 2 was relatively chaotic, while no strong flow field was formed at observation points 1 and 3. In addition, sound pressure levels of three observation points gradually decreased with the frequency increase.

Fig. 13Wind tunnel experiment of the rudimentary landing gear

Fig. 14Testing process of flow noise for the rudimentary landing gear

In order to verify accuracy of numerical simulation for flow noise of the rudimentary landing gear, a wind tunnel experiment was conducted, as shown in Fig. 13. The rudimentary landing gear was fixed in the wind tunnel, and an inflow with certain intensity was provided towards the rudimentary landing gear. Microphones were arranged at three points in Fig. 11. Pulse test system was used to record data. Each group of data was recorded for three times, and sampling frequency was 65536 Hz. The detailed test process was shown in Fig. 14. Sound pressure levels of the experimental were compared with the numerical simulation results, as shown in Fig. 15.

It can be found from Fig. 15 that sound pressure levels of the experimental values and that of the numerical simulation did not differ much, and the change trends were basically consistent. It was further shown in this result that the algorithm mentioned in the paper was feasible for computing the flow noise of the rudimentary landing gear.