Abstract
Sound field prediction has practical significance in the control of noise generated by sources in a flow, for example, the noise in aero-engines and ventilation systems. Aiming at accurate and flexible prediction of time-dependent sound field, a finite-difference wavenumber-time domain method for sound field prediction in a uniformly moving medium is proposed. The method is based on the second-order convective wave equation, and the wavenumber-time domain representation of the sound pressure field on one plane is forward propagated via a derived recursive expression. In this paper, the recursive expression is first deduced, and then numerical stability and dispersion of the proposed method are analyzed, based on which the stability condition is given and the correction of dispersion related to the transition frequency is made. Numerical simulations are conducted to test the performance of the proposed method, and the results show that the method is valid and robust at different Mach numbers.
1. Introduction
The noise generated by sources located in a moving medium, such as wind turbines, aero-engine rotors, and aircraft landing gears and airfoils, raises more and more concerns in recent years because such noise gives rise to negative effects on human beings and mechanical structures. In order to control this kind of noise, it is important to understand its propagation properties in the moving medium. Sound field prediction is exactly an effective way to comprehend such properties and thus can provide a basis for noise reduction.
Currently, several methods, such as parabolic equation approximations [1, 2], fast field program [3, 4], nearfield acoustic holography [5–7], and finite-difference time-domain (FDTD) method [8–10], have been presented for sound field prediction in a moving medium. The first two methods are performed in the frequency domain, while the rest can be performed in the time domain. Compared with the frequency domain method, the time domain method has some advantages that the transient effects can be easily observed, and the response over a broad frequency range can be obtained with a single simulation run on condition that the excitation is given as a short acoustic pulse [11, 12]. As for nearfield acoustic holography, it needs a known impulse response or green function in a moving medium, which can be specifically called real-time nearfield acoustic holography or time-domain equivalent source method, and in fact there is a distortion problem for impulse response sampling or an instability problem in iterative process to obtain equivalent source strength, respectively [7, 13]. On the contrast, the FDTD method does not need the transfer relationships, and thus it has advantages of flexibility and conceptual simplicity. Nevertheless, in the FDTD method, the wave propagation equations are directly solved in the space-time domain and the acoustic quantities at each spatial point are calculated in terms of those at nearby points. Therefore, the FDTD method is usually regarded as a local method [14]. Comparatively, some scholars have investigated the spectral methods (such as the k-space method and the pseudospectral method) [14–17], where the spatial Fourier transforms are used to evaluate the spatial derivatives in the wave propagation equations, and these methods are viewed as global because information from the entire wave field is employed to obtain the acoustic quantities at each spatial point. With the aid of spatial Fourier transforms, which can be implemented by fast Fourier transform algorithm, the global methods have advantages of greater accuracy and computational efficiency than the FDTD method. However, the developed spectral methods are based on the wave propagation equations in a stationary medium, and thus they cannot be used to realize sound field prediction in a moving medium, which is affected by the flow effects.
To realize the prediction of time-dependent sound field in a uniformly moving medium via a global method without any transfer relationship needed, a finite-difference wavenumber-time domain (FDWTD) method is proposed. The method is similar to but not the same as the spectral method. In the present method, the second-order convective wave equation is first transformed into the wavenumber-time domain by using the 2D spatial Fourier transform and then numerically solved by a finite difference scheme to obtain a recursive expression. By virtue of the recursive expression, if the time-dependent sound pressure on a plane in the nearfield of the sources is obtained, the pressure on any plane of interest further from the initial plane can be predicted, and thus the prediction of the three-dimensional sound field also can be realized easily. Due to the use of spatial Fourier transforms, the described method can be viewed as a global method. Furthermore, the method is derived from the convective equation; therefore, the flow effects caused by the moving medium are considered.
This paper is organized as follows. The basic theory of the proposed FDWTD method, the numerical stability, and dispersion analysis are presented in Section 2. In Section 3, numerical simulations are conducted to test the performance of the proposed method. Finally, conclusions are summarized in Section 4.
2. Theory
2.1. The FDWTD Method in a Uniformly Moving Medium
The second-order convective wave equation in a uniformly moving medium can be represented as [18]where , , and represent the sound velocity in a stationary medium, the emission time, and the sound pressure, respectively; represents the velocity of a moving medium; denotes the Laplace operator defined as ; is defined as , in which , , and represent the unit vectors in , , and directions, respectively.
Assume that the velocity of a moving medium is in the positive direction; then, (1) can be written aswhere is the magnitude of .
Taking the 2D spatial Fourier transform with respect to and , (2) can be expressed aswhere and represent the wavenumbers in and directions, respectively; represents the imaginary unit; represents the complex sound pressure in the wavenumber-time domain (here note that and for annotation in the bracket are omitted in view of brevity).
Equation (3) is the second-order convective wave equation in the wavenumber-time domain. Comparing (3) with (2), it can be found that the spatial derivatives with respect to and in (2) are evaluated using the spatial Fourier transform in (3). Unlike the finite difference stencils used in the FDTD method, in the present method the spatial derivatives with respect to and are based on the information of all points on the (x, y) plane instead of nearby points, and thus the present method is relatively global and yields a very high order spatial derivative approximation, equal to the number of grid points in each direction. This property has been proven by Fornberg [16].
Based on the central finite difference approximation, the first-order and second-order derivatives of the sound pressure with respect to time and displacement in the direction can be represented aswhere and represent the increments of time and grid interval in the direction, respectively, and represents the high order quantity that is often neglected in the finite difference approximation.
The substitution of (4) into (3) yieldswhere is the grid ratio defined as , is defined as with , and represents the Mach number defined as . Here, it is noted that when the medium flows along the positive direction, is positive, and negative conversely.
For the sake of simplicity, (5) can be rewritten aswhere and are the positive integers representing the numbers of spatial and time steps, respectively. Equation (6) is the recursive expression for calculating the wavenumber-time domain sound pressure on the plane and at the time from those on the adjacent planes and at the two previous time steps. The grid for calculation can be visually illustrated in Figure 1. To start the recursive calculation, it is given that represents the sound pressure on the initial plane to be propagated, and and are both equal to zero based on the assumption that there is no signal in the predicted sound field before and when the first signal sample is to be propagated. It is worthwhile to note that the sound pressure here is represented in the wavenumber-time domain, and (6) is the evolution of sound pressure at each wavenumber point .
The geometry of the proposed method is shown in Figure 2, and the procedure can be depicted as follows. Firstly, the wavenumber spectrum of the time-dependent sound pressure on the initial plane located at zh is calculated by 2D spatial Fourier transform. In practice, the initial plane is located at the nearfield of the sources, and the sound pressure on the plane could be measured or simulated. Then, by means of (6), the wavenumber-time domain sound pressure on the prediction plane located at zr can be obtained. Finally, by applying the inverse 2D spatial Fourier transform, the space-time domain sound pressure on the prediction plane is able to be acquired. In fact, the prediction plane can be any plane of interest further from the initial plane and consequently the three-dimensional sound field in a moving medium can be predicted.
2.2. Stability Analysis
To avoid numerical instability, the recursive computations require determinations of spatial and temporal sampling criteria. The stability of the scheme of (6) is analyzed based on the Von Neumann analysis [19]. According to references [20, 21], the analysis process can be implemented by employing the discrete Fourier transform and Z-transform. The wavenumber spectrum of pressure can be obtained from the spatial discrete Fourier transform:where is the wavenumber in direction.
It can be obtained from (7) that the wavenumber spectrum of pressure shifted by the spatial steps is able to be represented as
The substitution of (7) and (8) into (6) yieldswhere .
Taking the Z-transform with respect to time in (9), it can be represented asin which denotes the variable after Z-transform, and thus the characteristic or amplification equation can be obtained as
The stability condition of (6) is that, for all values of that is involved in the definition of a, the absolute roots of the amplification equation are less than or equal to unity. From (11), the inequality can be obtained as
To obtain the explicit expression of the stability condition, inequality (12) can be discussed in two different cases, according to the relationship between the magnitudes of and . It can be analyzed briefly as follows.
Case 1.
When the sign in inequality (12) takes “+,” the expression can be obtained asIt can be obviously found that inequality (13) could not be satisfied for all values of .
When the sign in inequality (12) takes “–,” it yieldswhich cannot be satisfied for each wavenumber in the case of taking equality in inequality (14).
Case 2.
Under this condition, whether the sign in inequality (12) takes “+” or “–,” it yields the same expression as follows:This inequality can be always satisfied unconditionally.
Based on the analysis above, it can be concluded that can be unconditionally satisfied when , while cannot be always satisfied when . Thus, the expression of stability condition constrained by inequality (12) isSubstituting into inequality (16) and applying the trigonometric formula, it yieldsThis condition should be satisfied for all values of , and thus it can be rewritten asSubstituting into inequality (18), the stability condition eventually can be represented aswhere is defined as and denotes the sampling frequency.
2.3. Dispersion Analysis
Different from the stability analysis that is conducted in the wavenumber domain, the dispersion analysis should be performed in both the frequency and wavenumber domains because the dispersion is inherently related to frequency. The frequency spectrum of can be obtained aswhere represents the angular frequency.
It can be obtained from (20) that the frequency spectrum of shifted by the time steps is able to be represented as
Substituting (20) and (21) into (9), it yields
Here it is worthy to note that, from the definition following (9), is an even function of the wavenumber , and thus (22) should be actually represented aswhere the sign “” represents the absolute value and is defined as .
It has shown that should satisfy the condition in the stability analysis. Now it is important to choose a proper specific value for , which will determine the dispersion relationship and further determine the transition frequency that separates two kinds of behavior of acoustic waves: propagating waves and evanescent waves. The transition frequency will be deviated from that in the continuous case caused by discretization, and it is expected to be corrected in the following analysis.
Firstly assuming that and according to inequality (18), the relationship between and can be obtained as
Combining (23) and (24) and the definition of , the dispersion relationship between the angular frequency and the wavenumber can be obtained as
Substituting and into (25), the wavenumber in a uniformly moving medium normalized by can be yielded asin which is the acoustic wavenumber in the stationary medium defined as .
It can be found that there exists a transition frequency at which changes from the imaginary value to the real value, and the frequency normalized by the sampling frequency can be represented as
In fact, the theoretical value of wavenumber in the continuous case normalized by in the propagation direction can be represented aswhere [5]. From (28), the normalized transition frequency can be obtained as
From (27) and (29), it can be found that the transition frequency in the discrete case is deviated from that in the continuous case, and both of them depend on the velocity of the moving medium, the wavenumber point , and the sampling frequency. The deviation indicates that, in the normalized frequency band between and , the wavenumber in the propagation direction may become a real value due to the discretization while it should be an imaginary value actually. That is to say, at those frequencies the evanescent waves in the continuous case become the propagating waves in the discrete case. To illustrate the deviation more clearly, the comparisons between (26) and (28) at two different wavenumber points are shown in Figure 3. Here the flow velocity represented as the Mach number is given as , and the sampling frequency is chosen as 12.8 kHz. From Figure 3, the deviation of transition frequency can be observed, and it is more apparent at a higher wavenumber point. It is expected that the correction can be made by assuming that equals , and thus the corrected value of satisfying the assumption can be obtained as
(a)
(b)
(c)
(d)
The comparison at the wavenumber point (40, 40) after this correction is shown in Figure 4. It can be observed that the transition frequency in the discrete case has been shifted to the theoretical one.
(a)
(b)
So far, the case where has been considered. If , the relationship between and becomes , and a corresponding analysis can be conducted similarly to the previous process. However, compared with the case in which , a degradation of the results for can be observed. For instance, when and are chosen, the results at the wavenumber point (40, 40) before and after corrections are presented in Figure 5.
(a)
(b)
(c)
(d)
In fact, the dispersion relationship indicates that the dispersion is not only related to the transition frequency, but also related to the deviation of the magnitude of . However, it is difficult to make a correction to the magnitude because it is frequency dependent. Fortunately, after the correction of the transition frequency, the magnitude deviation is mainly found at high wavenumber and at high frequency, which is close to one half of the sampling frequency. Thus, it is possible to get around this problem by selecting a proper sampling frequency according to the analyzed frequency, which makes it a little more restricted than the Nyquist sampling theorem.
Based on the above analysis, the best configuration corresponds to . Substituting this condition and replacing with the corrected value , the recursive expression of the wavenumber-time domain method used from now on can be finally obtained as
3. Numerical Simulations
When the FDWTD method is applied for sound field prediction in an unbounded field, the absorbing boundary condition should be adopted to avoid the reflection because the computational region is limited in practice. However, it is enough here to choose a sufficiently large linear space so that it can satisfy the condition , where and represent the extremities of the analyzed region and observation time, respectively. Furthermore, the numerical stability and dispersion conditions require that , and thus the grid interval can be determined as
From (32), it is noted that changes with the values of and , which eventually depends on the wavenumber point . Therefore, if the propagation distance is given as , the number of steps that is necessary to reach the propagation distance for a certain wavenumber point can be determined by .
3.1. Parameters in the Simulations
Numerical simulations are carried out to verify the validity of the proposed method. The sound source consists of two monopoles located at S1 (0.1 m, 0.1 m, 0 m) and S2 (−0.1 m, −0.1 m, 0 m). The 2 m × 2 m initial plane with the grid interval of 0.05 m is located at 0.05 m. Two prediction planes are located at 0.2 m and 0.5 m. The sound pressure at can be calculated as [18]where and are represented aswhere is the coordinate of sound source, represents the derivative of with respect to time, and represents the strength of the source given by
The sampling frequency is set to be 51.2 kHz, and the length of the signal is 512 . To reduce the wraparound error associated with the 2D spatial Fourier transform, the sound pressure on the initial plane is zero padded, and the size after zero-padding is four times the previous one.
3.2. Results and Discussions
In the simulations, the velocity of the moving medium presented as the Mach number is first set to be 0.6. The sound pressures calculated by the proposed method on the above-mentioned two prediction planes at different time steps are compared with their theoretical values. Here, two time steps, 150 and 250 , are chosen. On each prediction plane, two points are selected to compare the calculated time signals with the theoretical ones, locating at P1 (0 m, 0 m, zr) and P2 (0.3 m, 0.8 m, zr). The theoretical and calculated pressures on the two prediction planes at the two time steps are shown in Figures 6 and 7. The calculated time signals at the two points on each plane are compared with the theoretical ones, shown in Figure 8.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
It can be seen from Figures 6 and 7 that the calculated pressures agree very well with the theoretical ones on each plane at the same time steps. By comparing the calculated time signals with the theoretical ones at the two points on each prediction plane, it can be observed from Figure 8 that they are also consistent with each other well, although a slight deviation occurs in the final time range. Furthermore, it can be found that the results at the central point P1 are slightly better than those at the point P2. These errors might be caused by the truncation of the sound field due to the limited aperture on the initial plane. In fact, the better performance can be observed if a larger initial plane is chosen in the simulations. Considering the brevity of the content and the errors being acceptable, the results for a larger initial plane are not given here. In practice, the size of the initial plane is a compromise of accuracy requirement and computational efficiency.
In order to evaluate the results more accurately and comprehensively, two temporal indicators T1 and T2 are defined bywhere denotes the time averaged value, is the calculated pressure by the proposed method, and is the theoretical pressure. T1 is a correlation coefficient which is sensitive to the similarity between the shapes of the signals and thus indicates their phase difference. The best value for T1 would be 1. T2 is the relative difference between the root mean square values of the two signals for characterizing the similarity of amplitudes. The target value to reach for T2 is 0. The spatial maps of the indicators are shown in Figure 9. From the results, it can be seen that the values of indicator T1 are very close to 1 in both cases, which means that the phases of the calculated pressure and theoretical one are very similar to each other, and the values of indicator T2 are small (less than 0.2 overall in both cases), which means that the differences between the magnitudes of the calculated pressure and theoretical one are small. Although some deviations can be observed in the corner or edge areas, which might be caused by the limited aperture as mentioned above, the results indicate that the calculated pressure is satisfactory compared with the theoretical one in both spatial meaning and temporal meaning.
(a)
(b)
(c)
(d)
Furthermore, the indicators are also calculated at other Mach numbers on plane 0.2 m, to check the performances of the method at different flow velocities. In order to make the results easily compared, the arrays of the indicator values are rearranged into lines in the order from up to down and from left to right, and the results are shown in Figure 10. It can be found that the values of indicator T1 at different Mach numbers are close to 1, more than 0.9 for all points. Meanwhile, the values of indicator T2 at different Mach numbers are at a low level, less than 0.2 except one or two points. In addition, it also can be seen that there are a few points in worse performance than others with a certain period, and this phenomenon is caused by the permutation of the indicator values and the error due to limited aperture mentioned above. The results indicate that the prediction errors are not amplified by the increasing flow effects as the velocity of the moving medium increases. This is because the flow effects due to the moving medium have been fully considered in the proposed method, which can be reflected in the second-order convective wave equation, the stability condition, and the correction of dispersion related to transition frequency.
(a)
(b)
4. Conclusions
A FDWTD method for sound field prediction in a uniformly moving medium has been proposed based on the second-order convective wave equation. In the method, the second-order convective wave equation is first transformed into the wavenumber-time domain by taking the 2D spatial Fourier transform, then the first-order and second-order derivatives of sound pressure in the wavenumber-time domain are approximated by the central finite difference in second order, and finally a recursive expression is deduced for the sound field prediction. Considering the numerical stability and dispersion in the calculation procedure, the stability condition has been given with a full consideration of the moving medium, and the correction of dispersion related to the transition frequency has been made based on the assumption that the transition frequency in the discrete case should be equal to that in the continuous case.
Numerical simulations have been conducted to verify the validity and robustness of the proposed method, and the results have shown that the calculated sound pressures are in good agreement with the theoretical ones and the proposed method works stably at different Mach numbers. Furthermore, it should be noted that the described method has the potential for sound field prediction in a stratified moving medium. In addition, the absorbing boundary condition, as a widely studied area [22–24], is worthwhile to be extended to address the problem of boundary effects in the present method, and it is expected to be studied in the future.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant no. 11804002), the University Science Research Project of Anhui Province (grant no. KJ2019A0792), and the Anhui Jianzhu University Research Project (grant no. 2018QD06).