Finite-Difference Frequency-Domain Scheme for Sound Scattering by a Vortex with Perfectly Matched Layers

Abstract

Understanding the effect of vortexes on sound propagation is of great significance in the field of target detection and acoustic imaging. A prediction algorithm of the two-dimensional vortex scattering is realized based on a finite-difference frequency-domain (FDFD) numerical scheme with perfectly matched layers (PML). Firstly, the governing equation for flow–sound interaction is given based on the perturbation theory, and the FDFD program is built. Subsequently, the mesh independence is verified, and the result has a good convergence when the mesh corresponds to over 15 nodes per wavelength. Then, computational parameters of the PML are discussed to achieve better absorbing boundary conditions. Finally, the results of this algorithm are compared with previous literature data. Results show that for different cortex scattering cases, the absorption coefficient should vary linearly with the density of the medium and the incident wave frequency. When the thickness of the PML boundary is greater than 2.5 times the wavelength, the PML boundary can absorb the scattering sound effectively. This provides a reliable algorithm for the numerical study of the effect of vortexes on sound propagation.

1. Introduction

When sound waves propagate through vortical flow, the flow–sound interaction affects the sound propagation [1], resulting in the acoustic scattering phenomena [2,3,4]. Two-dimensional vortex scattering is a classical model in the understanding of flow–sound interaction [5,6,7,8,9]. The study of the characteristics of vortex scattering is important for both exploring the coupling mechanism of flow and acoustic [10] and capturing the features of the sound propagation in nonlinear flow in engineering, such as sound propagation through the wake of the aircraft jets and acoustic imaging in medicine.

Numerical algorithms are becoming more accurate [11,12,13,14,15]. Numerical simulation methods are the main approach to deal with the problem of vortex acoustic scattering [16,17]. Coloninus [18] analyzed the propagation of acoustic waves through a Gaussian vortex flow based on a high-order numerical scheme. Karabasov [19] solved the Euler equation with the extended upwind leapfrog scheme and analyzed the acoustic scattering from steady-state Gaussian vortex and Rankine vortex flow at different Mach numbers. Blanc-Benon [20] used geometric acoustic theory and parabolic wave equations to simulate the propagation of acoustic waves through the turbulent flow. Iwatsu and Tsuru [21] used Cin’s integral method and the compact finite-difference scheme to solve the linear Euler equations and calculated the scattering acoustic field from the Rankine vortex and Burgers vortex, respectively. Ke [2] analyzed the acoustic scattering characteristics of plane waves propagating through the Gaussian vortex and vortex pairs with the WENO scheme, and the results were consistent with the direct numerical simulation obtained by Coloninus [18]. In the study of underwater vortex scattering, the finite-difference time-domain method (FDTD) was applied [22,23,24]. Zhang [25] introduced the simulation of vortex scattering based on the FDTD method and discussed the directivity of the scattering sound.

To improve the stability in the numerical simulation of vortex scattering, the finite-difference frequency-domain (FDFD) algorithm with numerical discretization on the staggered grid is introduced to the simulation of vortex scattering.

Meanwhile, the construction of an acoustic absorption boundary is especially important for numerical simulation. By adding several absorption layers around the computational domain, the perfectly matched layers (PML) can decrease the boundary reflection and thus turn the infinite domain into a finite domain. For acoustic propagation in the moving flow, PML can be used as the acoustic absorbing boundary [26]. The acoustic absorption effect of PML depends on the absorption coefficient, and setting a proper absorption coefficient improves the acoustic absorption [27]. Moreira [28] solved the propagation equation of acoustic waves in uniform and inhomogeneous medium using the frequency-domain finite-difference method combined with PML and analyzed the influence of the PML absorption coefficient on the acoustic absorption effect. Jie [29] adopted PML for the broadband scattering model of the underwater target and set different absorption coefficients according to the frequency of the incident wave. It shows that the acoustic absorption effect works well at each frequency. Park [30] applied the PML to solve the 2D Helmholtz equation and symmetrically discretized the PML matrix, which can effectively reduce the computational resource. However, the effect of PML on the calculation of vortex scattering and the suitable value of parameters for vortex scattering simulation are unclear, such as the effect of the relative thickness of PML, the size of the computational domain, and the vortex core radius. In order to study the above issues quantitatively, we have introduced an attenuation term in the equation of motion in this paper, which increases from the inner to the outer boundary, instead of the diffusion term multiplied by the attenuation factor according to the aforementioned works. This method can ensure energy absorption, and the attenuation variation inside the PML layers is relatively smooth. This ensures the precision of numerical simulations regarding the scattering of vortex acoustics.

In this paper, a FDFD algorithm combined with the PML is established to evaluate the acoustic scattering caused by a vortex, and the relevant parameters suitable for 2D vortex acoustic scattering are specially discussed and determined. In the end, the validated results agree well with the reference results. In Section 1, a two-dimensional vortex scattering model is introduced. In Section 2, the numerical scheme to solve the perturbation equation is given. In Section 3, the algorithm is introduced, and the grid independence is verified. In Section 4, the relationship between the PML absorption coefficient, PML layers, and Mach number is explored, and the algorithm is verified. Section 5 contains the conclusion.

2. Vortex Scattering Model

In this paper, the classical two-dimensional Burgers vortex model is used. The schematic diagram of vortex scattering is shown in Figure 1:

The tangential velocity of the vortex flow is [2]:

$$v_{\theta }=\frac{\Gamma }{2\pi }\left[1-\mathit{exp}(-\beta r^2/a^2)\right],$$

(1)

where $a$ is the radius of the vortex core. $r$ is the distance from the point in the flow to the center of the vortex core. $\beta$ = 1.256431 is a constant chosen so that the maximum velocity occurs at $r/a=1$. $\Gamma =2.8\pi L{c}_{0}Ma$ is the circulation of the vortex. The Mach number $Ma={(v_{\theta })}_{\max }/c_0$ is the ratio of the maximum velocity to the sound speed.

The sound pressure of the incident plane wave is defined as [31]:

$$p_{in}=p_a\mathit{sin}(2\pi ft-kx),$$

(2)

where $p_a$ is the amplitude of the incident wave, $f$ is the frequency, $k$ is the wavenumber, $x$ is the $x$-direction coordinate, and t is the time.

3. Governing Equations

For vortex flow with a low Mach number, the amplitude of acoustic variables is assumed to be smaller than the fluid variables. Therefore, based on the first-order perturbation method, the Navier–Stokes equations can be separated into two parts based on the magnitude. Then, the wave equation of the flow–sound interaction is obtained.

For flow without external force, according to the isentropic adiabatic principle, the governing equations can be simplified [32]:

$$\frac{\partial p}{\partial t}=-\rho c^2\nabla \cdot \mathit{u}-(\mathit{u}\cdot \nabla )p,$$

(3)

$$\frac{\partial \mathit{u}}{\partial t}=-\frac{1}{\rho }\nabla p+\upsilon {\nabla }^2\mathit{u}-(\mathit{u}\cdot \nabla )\mathit{u},$$

(4)

where $\rho$, $\mathit{u}$, and $p$ represent the density, velocity and pressure of the flow, respectively. $\upsilon$ is the kinematic viscosity, and c is the sound speed.

Assuming that the magnitude of the acoustic variable is much smaller than the corresponding fluid variable [33], p and u can be defined as follows:

$$\{\begin{array}{l}p=p_0+\delta p,\left|\delta p\right|<<p_0\\ \mathit{u}={\mathit{u}}_0+\delta \mathit{u},\left|\delta \mathit{u}\right|<<\left|{\mathit{u}}_0\right|\end{array},$$

(5)

where $p_0$ and ${\mathit{u}}_{\mathbf{0}}$ represent the pressure and velocity of the flow without the acoustic disturbance, respectively. $\delta p$ and $\delta \mathit{u}$ are the perturbed acoustic pressure and velocity generated by the sound wave.

We consider $\delta p=\delta \rho c^2$, so that $p_0$ and ${\mathit{u}}_{\mathbf{0}}$ satisfy the following formulas:

$$\frac{\partial p_0}{\partial t}=-{\rho }_0{c}^2\nabla \cdot {\mathit{u}}_0-({\mathit{u}}_0\cdot \nabla )p_0,$$

(6)

$$\frac{\partial {\mathit{u}}_0}{\partial t}=-\frac{1}{{\rho }_0}\nabla p_0+\upsilon {\nabla }^2{\mathit{u}}_0-({\mathit{u}}_0\cdot \nabla ){\mathit{u}}_0,$$

(7)

For low Mach number flow, the following assumptions are met: (1) $\left|\delta \mathit{u}\right|<<\left|{\mathit{u}}_0\right|<<c$, (2) the time scale of the flow is much larger than the period of the incident acoustic wave, and (3) the influence of the acoustic disturbance on the flow is negligible. The vortex scattering phenomenon caused by the flow–sound interaction is then mainly controlled by the nonlinear term. By retaining the term and neglecting the high-order amount, the simplified governing equation of acoustic disturbance can be obtained:

$$\frac{\partial \delta p}{\partial t}=-{\rho }_0{c}^2\nabla \cdot \delta \mathit{u}-({\mathit{u}}_0\cdot \nabla )\delta p,$$

(8)

$$\frac{\partial \delta \mathit{u}}{\partial t}=-\frac{1}{{\rho }_0}\nabla \delta p-({\mathit{u}}_0\cdot \nabla )\delta \mathit{u}-(\delta \mathit{u}\cdot \nabla ){\mathit{u}}_0,$$

(9)

where ${\rho }_0$ is the density of the medium.

The above time domain equation can be transformed into a frequency domain equation. It can be assumed that $\delta p$ and $\delta \mathit{u}$ satisfy $\delta p=\delta p_a\cdot e^{j\omega t}$, $\delta \mathit{u}=\delta {\mathit{u}}_a\cdot e^{j\omega t}$ respectively, and then Equations (8) and (9) can be transformed into:

$$j\omega \delta p_a=-{\rho }_0{c}^2\nabla \cdot \delta {\mathit{u}}_a-({\mathit{u}}_0\cdot \nabla )\delta p_a,$$

(10)

$$j\omega \delta {\mathit{u}}_a=-\frac{1}{{\rho }_0}\nabla \delta p_a-({\mathit{u}}_0\cdot \nabla )\delta {\mathit{u}}_a-(\delta {\mathit{u}}_a\cdot \nabla ){\mathit{u}}_0,$$

(11)

where $\omega$ is the circular frequency.

The PML is adopted as the acoustic absorption boundary, as shown in Figure 2. The computational domain is a square with N-layers PML. In order to effectively absorb the sound, the attenuation factor of the i-th layer satisfies the following formula:

$$\sigma (i)={\sigma }_{\max }{(\frac{i-1}{N-1})}^2,$$

(12)

where ${\sigma }_{\max }$ is the maximum absorption coefficient.

Obviously, in the PML, Equation (11) can be modified to:

$$j\omega \delta {\mathit{u}}_a=-\frac{1}{{\rho }_0}\nabla \delta p_a-({\mathit{u}}_0\cdot \nabla )\delta {\mathit{u}}_a-(\delta {\mathit{u}}_a\cdot \nabla ){\mathit{u}}_0-\frac{\sigma }{{\rho }_0}\delta {\mathit{u}}_a,$$

(13)

It can be seen from Equation (13) that the attenuation term adds an attenuation factor to the velocity of the particle inside the PML layers.

4. Numerical Method and Mesh Independence

4.1. Numerical Method for Acoustic Scattering

To ensure the accuracy of the simulation, the finite difference scheme with the staggered grid is applied to solve Equations (10) and (11). Sound pressure $\delta p$ corresponds to each staggered node. The velocity components in the $x$ and $y$ direction, respectively, correspond to the node with $∆x$/2 offset in the $x$-direction and $∆y$/2 offset in the y-direction from the pressure node. The two-dimensional staggered mesh is shown in Figure 3.

Equations (10) and (11) can be rewritten as:

$$\left\{\begin{array}{c}\frac{\partial {\delta u}_{ax}}{\partial x}=\frac{1}{dx}({\delta u}_{ax}(i,j)-{\delta u}_{ax}(i-1,j)),\\ \frac{\partial {\delta u}_{ay}}{\partial y}=\frac{1}{dy}({\delta u}_{ay}(i,j)-{\delta u}_{ay}(i,j-1)),\\ \frac{\partial \delta p_a}{\partial x}=\frac{1}{dx}(\delta p_a(i,j)-{\delta p}_a(i-1,j)),\\ \frac{\partial {\delta p}_a}{\partial y}=\frac{1}{dy}(\delta p_a(i,j)-{\delta p}_a(i,j-1)),\end{array}\right.$$

(14)

where ${\delta u}_{ax}$ and $\delta u_{ay}$ are the amplitudes of velocity in x- and y-axis direction, respectively, $\delta p_a$ is the amplitude of sound pressure. Moreover, $\frac{\partial {\delta p}_a}{\partial x}$ and $\frac{\partial \delta p_a}{\partial y}$ are the partial differential of sound pressure.

4.2. Mesh Independence

This section performs the study of mesh independence. We compare the results from different mesh with the different number of grid nodes (10, 15, and 20) per wavelength. The sound scattering directivity results of these cases are compared in Figure 4. ${\delta p}_{rms}$ is the effective scattered sound pressure, $\theta$ is the azimuth. It can be seen from the figure that the solution has a very good convergence when the mesh corresponds to over 15 nodes per wavelength.

5. Discussion of PML Parameters

5.1. Relationship between Layer Numbers and Absorption Coefficient

Firstly, the effect of the absorption coefficient on the result for a given number of PML layers is discussed. By choosing sufficient thickness of PML as 75 layers (5 times the wavelength) and changing ${\sigma }_{max}$ from 200 to 30,000 ${\mathrm{kg}/(\mathrm{m}}^3\cdot \mathrm{s})$, the sound scattering directivity results of different cases are compared in Figure 5.

We define the relative error as follows:

$$err=\frac{{p^{\prime }}_m-p_m}{p_m},$$

(15)

where ${p^{\prime }}_m$ is the effective scattered sound pressure at the monitoring point $m$ of the test. $p_m$ is the theoretical value of effective scattered sound pressure at the monitoring point $m$.

As can be seen from Figure 5, obvious fluctuations in the backward scattering region arise with an absorption coefficient too small or too large. When the attenuation is too small (200 kg/(m³·s), the PML boundary absorbs less acoustic energy, resulting in excessive reflection. When the attenuation is too large, the boundary behaves like a rigid boundary and results in a stronger reflection. Therefore, an appropriate absorption coefficient should be chosen.

In order to analyze the attenuation of sound pressure inside the PML layers, the distribution of effective scattering sound pressure of the PML boundary is shown in Figure 6. When ${\sigma }_{max}$ = 200 ${\mathrm{kg}/(\mathrm{m}}^3\cdot \mathrm{s})$, the scattered sound pressure is not completely attenuated in PML layers, and the distribution along the PML layers shows significant fluctuations. When ${\sigma }_{max}$ = 1000–15,000 ${\mathrm{kg}/(\mathrm{m}}^3\cdot \mathrm{s})$, the scattered acoustic pressure in PML layers can be attenuated smoothly, and the larger the absorption coefficient is, the fast it attenuates. When ${\sigma }_{max}$ = 30,000 ${\mathrm{kg}/(\mathrm{m}}^3\cdot \mathrm{s})$, the scattered acoustic pressure can be completely attenuated in PML layers, but local fluctuations occur near the first layer of the PML layers, which produces a strong reflection into the internal computational domain.

It can be seen that the appropriate absorption coefficient should satisfy the following conditions: the scattered acoustic pressure inside the PML layers can be completely attenuated, and no strong reflection is produced in the internal computational domain.

Another problem that should be considered is that too many PML layers bring low computing efficiency. Figure 7 shows the variation in the calculation time with the number of PML layers. As can be seen from the figure, when the number of PML layers increases, the calculation time also increases rapidly.

Based on this study of this section, when a larger number of layers is chosen, a better acoustic absorption can be obtained, but it may have a larger calculation load. Therefore, a reasonable number of layers should be chosen, which will be analyzed shortly.

5.2. Dimensionless Analysis of PML Parameters

This section analyzes the relationship between the absorption coefficient and related parameters based on the dimensionless analysis method. There are eight related physical quantities, namely medium density ${\rho }_0$, sound speed $c_0$, incident wave frequency $f$, number of nodes per wavelength $n$, the size of the computational domain $L$, Mach number $Ma$, vortex core radius $a$, number of PML layers $N$, and PML absorption coefficient ${\sigma }_{\max }$. We choose ${\rho }_0, f, c_0$ as the basic quantity. the dimensional relationship is:

$$\frac{{\sigma }_{\max }}{{\rho }_{0}f}=g(\frac{N}{n},\frac{Lf}{c_0},\frac{af}{c_0},Ma),$$

(16)

where $g(\frac{N}{n},\frac{Lf}{c_0},\frac{af}{c_0},Ma)$ is a function. $\frac{N}{n}$ is the ratio of the thickness of the PML boundary to the wavelength. $\frac{Lf}{c_0}$ is the ratio of the size of the computational domain to the wavelength. $\frac{af}{c_0}$ is the ratio of the radius of the vortex to the wavelength.

The absorption coefficient relationship is (17):

$${\sigma }_{\max }={\rho }_{0}fg(\frac{N}{n},\frac{Lf}{c_0},\frac{af}{c_0},Ma),$$

(17)

In Equation (17), $\frac{N}{n}$ is the main parameter characterizing the thickness of the PML boundary. From Equation (1), it can be seen that the flow velocity distribution is determined using $L,Ma,a$, therefore $\frac{Lf}{c_0}$, $\frac{af}{c_0}$, and $Ma$ can be regarded as the main parameters characterizing the flow velocity.

5.3. Parameters of Flow

In this section, the relationship between the absorption coefficient and the flow field velocity parameters $\frac{Lf}{c_0}$, $\frac{af}{c_0}$, and $Ma$ are analyzed. Since a large $L$ causes a large amount of calculations, we first analyze $\frac{Lf}{c_0}$, then $\frac{af}{c_0}$, and finally $Ma$ with proper $L$. The frequency is set to be 85 Hz, the sound speed is 340 m/s, the number of nodes per wavelength is 15, and ${\sigma }_{\max }=5000\mathrm{kg}/({\mathrm{m}}^3\cdot \mathrm{s})$.

5.3.1. Size of the Computational Domain

This section discusses the influence of the size of the computational domain on the simulation. The effective scattering attenuation distribution in the PML with different sizes of calculation domains (2.5, 5.0, 10.0) is shown in Figure 8. Results show that sound pressure is completely attenuated within the same number of layers, and there is no obvious local fluctuation. The size of the computational domain has little effect on the choice of PML parameters. Therefore, the following analysis of PML parameters can use a smaller calculation domain to reduce computational cost.

5.3.2. Vortex Core Radius

This section analyzes the radius of the vortex core. The incident sound frequency is 85 Hz, the sound speed is 340 m/s, the size of the computational domain is 10 m, and the radius of the vortex core is from 0.5 to 1.5 m. Therefore $\frac{af}{c_0}$ is from 0.125 to 0.375. The effective scattering sound distribution in PML layers is shown in Figure 9. It can be seen from the figure that the acoustic pressure scattered from different vortex is completely attenuated within the same number of layers, and there is no obvious local fluctuation. The vortex core radius has little effect on the choice of PML parameters.

5.3.3. Mach Number

This section analyzes the relationship between the Mach number and the absorption coefficient. The frequency is 85 Hz, the sound speed is 340 m/s, the size of the computational domain is 10 m, and the Mach number is 0.10–0.25. The scattering sound distribution in PML layers is shown in Figure 10. It can be seen from the figure that the acoustic pressure is completely attenuated in the same number of layers at different Mach numbers, and there is no obvious local fluctuation. The Mach number has little effect on the choosing of PML parameters.

It can be seen from the results in this section that the flow velocity parameters $\frac{Lf}{c_0}$, $\frac{af}{c_0}$, and $Ma$ have little effect on the parameter choosing of the PML boundary. When the flow velocity increases, the particle velocity increases. If the PML attenuation factor does not change, the scattered acoustic is attenuated at the same attenuation ratio in the PML layers. At the same time, due to the large attenuation factor, the PML absorption coefficient approximately satisfies the following formula:

$${\sigma }_{\max }={\rho }_{0}fg(\frac{N}{n}),$$

(18)

5.4. Relative Thickness of PML

This section analyzes the effect of PML thickness on acoustic absorption. The number of nodes per wavelength is chosen to be 15, and the number of PML layers $N$ is from 15 to 75. Therefore, $\frac{N}{n}$ is from 1.0 to 5.0. For different thicknesses, according to the function g ($\frac{N}{n}$), cases of different absorption coefficients are conducted. The comparison of the directivity of different cases is shown in Figure 11. It can be seen from the figure that when the thickness of PML is smaller than 2.5 times the wavelength, there is no obvious reflection.

In order to ensure effective acoustic absorption while minimizing the amount of calculation resource, a PML thickness of 2.5 times the wavelength is suggested.

5.5. Comparison with Direct Numerical Simulation (DNS)

In this section, the results obtained by our algorithm are compared with DNS results [18]. The density of the fluid medium is ${\rho }_0=1.0\hspace{0.33em}\mathrm{kg}/{\mathrm{m}}^3$, the characteristic radius of the vortex core is $L=1.0\mathrm{m}$, the Mach number is $Ma=0.25$, the frequency of the incident plane wave is $f=85\mathrm{Hz}$, the sound speed is $c=340\mathrm{m}/\mathrm{s}$, and the computational domain is a $40\mathrm{m}\times 40\mathrm{m}$ square. The monitoring points are $10a$ away from the center of the vortex core.

From Section 5.4, we know that when the PML boundary thickness is 2.5 times the wavelength, the reasonable $g(\frac{N}{n})$ is 23.5–40. The scattering sound pressure directivity compared with the DNS result is shown in Figure 12. It can be seen from the figure that the results of this method are consistent with the results obtained through DNS [18]; thus, the algorithm can be used to predict the vortex scattering field.

6. Conclusions

In this paper, a finite-difference algorithm in the frequency domain combined with the PML absorbing boundary is developed to simulate acoustic scattering from a two-dimensional vortex. The mesh independence is verified, and when the mesh corresponds to over 15 nodes per wavelength, the result has a good convergence. The value of the PML absorption coefficient and the number of PML layers are analyzed. The limitation of this algorithm is that it is only applicable to 2D vortex acoustic scattering. But it also provides a reliable algorithm for the numerical study of the effect of the vortex on sound propagation, such as target detection and acoustic imaging.

By comparing with the literature data, the accuracy of the algorithm is verified, and the following conclusions are drawn:

• In the PML region, when the scattered sound pressure can be completely attenuated before reaching the outermost boundary and there is no obvious fluctuation near the boundary between the computational domain and the PML, the scattering sound can be regarded as fully absorbed.
• When the thickness of the PML boundary is larger than 2.5 times the wavelength, a qualified sound absorption effect can be obtained by selecting a suitable PML absorption coefficient.
• The algorithm established in this paper can effectively calculate the two-dimensional vortex scattering. The value of the absorption coefficient of PML is suggested to be linear with the density of the medium and the frequency of the incident wave. The Mach number shows little effect on the acoustic absorption effect of PML.