A Finite-Difference Wavenumber-Time Domain Method for Sound Field Prediction in a Uniformly Moving Medium

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 domainmethod for sound field prediction in a uniformlymovingmedium is proposed.,emethod 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.


Introduction
e 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][6][7], and finite-difference time-domain (FDTD) method [8][9][10], have been presented for sound field prediction in a moving medium. e 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. erefore, 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][15][16][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. e 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.
is paper is organized as follows. e 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.

e FDWTD Method in a Uniformly Moving Medium.
e second-order convective wave equation in a uniformly moving medium can be represented as [18] 1 where c, t, and p represent the sound velocity in a stationary medium, the emission time, and the sound pressure, respectively; V represents the velocity of a moving medium; ∇ 2 denotes the Laplace operator defined as ∇ 2 � z 2 /zx 2 + z 2 /zy 2 + z 2 /zz 2 ; ∇ is defined as ∇ � z/zxi + z/zyj + z/zzk, in which i, j, and k represent the unit vectors in x, y, and z directions, respectively.
Assume that the velocity of a moving medium is in the positive x direction; then, (1) can be written as where V is the magnitude of V.
Taking the 2D spatial Fourier transform with respect to x and y, (2) can be expressed as where k x and k y represent the wavenumbers in x and y directions, respectively; i represents the imaginary unit; P represents the complex sound pressure in the wavenumbertime domain (here note that k x and k y 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 x and y 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 x and y 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. is 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 z direction can be represented as where Δt and Δz represent the increments of time and grid interval in the z direction, respectively, and O(·) represents the high order quantity that is often neglected in the finite difference approximation.
where r is the grid ratio defined as r � cΔt/Δz, b is defined x + k 2 y , and M represents the Mach number defined as M � V/c. Here, it is noted that when the medium flows along the positive x direction, V is positive, and negative conversely.
For the sake of simplicity, (5) can be rewritten as where m and n are the positive integers representing the numbers of spatial and time steps, respectively. Equation (6) is the recursive expression for calculating the wavenumbertime domain sound pressure on the plane mΔz and at the time (n + 1)Δt from those on the adjacent planes and at the two previous time steps. e grid for calculation can be visually illustrated in Figure 1. To start the recursive calculation, it is given that P(0, n) represents the sound pressure on the initial plane to be propagated, and P(m, 0) and P(m, 1) 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 P(0, 1) 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 (k x , k y ). e 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 z h 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. en, by means of (6), the wavenumber-time domain sound pressure on the prediction plane located at z r 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.

Stability Analysis.
To avoid numerical instability, the recursive computations require determinations of spatial and temporal sampling criteria. e 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. e wavenumber spectrum P(k z , n) of pressure P(m, n) can be obtained from the spatial discrete Fourier transform: where k z is the wavenumber in z direction. It can be obtained from (7) that the wavenumber spectrum P ±l (k z , n) of pressure P(m ± l, n) shifted by the spatial steps ±l is able to be represented as P ±l k z , n � P k z , n e ±ik z lΔz . (8) e substitution of (7) and (8) into (6) yields where a � 1 − r 2 − 2b 2 + r 2 cos k z Δz.
Known pressure Target pressure Sound field

Shock and Vibration
Taking the Z-transform with respect to time in (9), it can be represented as in which Z denotes the variable after Z-transform, and thus the characteristic or amplification equation can be obtained as e stability condition of (6) is that, for all values of k z 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 a 2 and 1 + V 2 k 2 x Δt 2 . It can be analyzed briefly as follows.
x Δt 2 When the sign in inequality (12) takes "+," the expression can be obtained as It can be obviously found that inequality (13) could not be satisfied for all values of k z .
When the sign in inequality (12) takes "-," it yields which cannot be satisfied for each wavenumber k z in the case of taking equality in inequality (14).
x Δt 2 Under this condition, whether the sign in inequality (12) takes "+" or "-," it yields the same expression as follows: is inequality can be always satisfied unconditionally. Based on the analysis above, it can be concluded that |Z 1,2 | ≤ 1 can be unconditionally satisfied when x Δt 2 , while |Z 1,2 | ≤ 1 cannot be always satisfied when a 2 > 1 + V 2 k 2 x Δt 2 . us, the expression of stability condition constrained by inequality (12) is Substituting a � 1 − r 2 − 2b 2 + r 2 cos k z Δz into inequality (16) and applying the trigonometric formula, it yields is condition should be satisfied for all values of k z , and thus it can be rewritten as (18), the stability condition eventually can be represented as where F t is defined as F t � cK M /2π and f s denotes the sampling frequency.

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. e frequency spectrum of P(k z , n) can be obtained as where ω represents the angular frequency. It can be obtained from (20) that the frequency spectrum P ±g (k z , ω) of P(k z , n ± g) shifted by the time steps ±g 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), a is an even function of the wavenumber k x , and thus (22) should be actually represented as where the sign "| · |" represents the absolute value and φ is defined as φ � arctan(1/V|k x |Δt). It has shown that a 2 should satisfy the condition a 2 ≤ 1 + V 2 k 2 x Δt 2 in the stability analysis. Now it is important to choose a proper specific value for a 2 , 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. e 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 a 2 � 1 + V 2 k 2 x Δt 2 and according to inequality (18), the relationship between b and r can be obtained as Combining (23) and (24) and the definition of a, the dispersion relationship between the angular frequency ω and the wavenumber k z can be obtained as Substituting b � πF t /f s and Δz � cΔt/r into (25), the wavenumber k z in a uniformly moving medium normalized by k 0 can be yielded as in which k 0 is the acoustic wavenumber in the stationary medium defined as k 0 � ω/c � 2πf/c. It can be found that there exists a transition frequency f t at which k z /k 0 changes from the imaginary value to the real value, and the frequency f t normalized by the sampling frequency f s can be represented as In fact, the theoretical value of wavenumber k z in the continuous case normalized by k 0 in the propagation direction can be represented as . From (28), the normalized transition frequency f t /f s can be obtained as (29) 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 (k x , k y ), and the sampling frequency. e deviation indicates that, in the normalized frequency band between f t /f s and f t /f s , the wavenumber in the propagation direction may become a real value due to the discretization while it should be an imaginary value actually. at 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 M � 0.6, 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 f t /f s equals f t /f s , and thus Shock and Vibration 5 the corrected value of F t /f s satisfying the assumption can be obtained as e comparison at the wavenumber point (k x , k y ) � (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.
So far, the case where a 2 � 1 + V 2 k 2 x Δt 2 has been considered. If a 2 < 1 + V 2 k 2 x Δt 2 , the relationship between b and r becomes b 2 + r 2 < 1/2(1 + ���������� 1 + V 2 k 2 x Δt 2 ), and a corresponding analysis can be conducted similarly to the previous process. However, compared with the case in which a 2 � 1 + V 2 k 2 x Δt 2 , a degradation of the results for k z /k 0 can be observed. For instance, when a 2 � 1/4 and b 2 � 3/4 − r 2 are chosen, the results at the wavenumber point (k x , k y ) � (40, 40) before and after corrections are presented in Figure 5.
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 k z /k 0 . 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. us, 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 b with the corrected value b � πF t /f s , the recursive expression of the wavenumber-time domain method used from now on can be finally obtained as

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 z max and T 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 Δz can be determined as From (32), it is noted that Δz changes with the values of b and k x , which eventually depends on the wavenumber point (k x , k y ). erefore, if the propagation distance is given as d, the number of steps N k that is necessary to reach the propagation distance for a certain wavenumber point can be determined by N k � d/Δz.

Parameters in the Simulations.
Numerical simulations are carried out to verify the validity of the proposed method. e sound source consists of two monopoles located at S 1 (0.1 m, 0.1 m, 0 m) and S 2 (− 0.1 m, − 0.1 m, 0 m). e 2 m × 2 m initial plane with the grid interval of 0.05 m is located at z h �0.05 m. Two prediction planes are located at z r �0.2 m and 0.5 m. e sound pressure at r � (x, y, z) can be calculated as [18] where R and R 1 are represented as where (x 0 , y 0 , z 0 ) is the coordinate of sound source, q ′ (t) represents the derivative of q(t) with respect to time, and q(t) represents the strength of the source given by e sampling frequency is set to be 51.2 kHz, and the length of the signal is 512 Δt. 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.

Results and Discussions.
In the simulations, the velocity of the moving medium presented as the Mach number is first set to be M �0.6. e 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 Δt and 250 Δt, are chosen. On each prediction plane, two points are selected to compare the calculated time signals with the theoretical ones, locating at P 1 (0 m, 0 m, z r ) and P 2 (0.3 m, 0.8 m, z r ). e theoretical and calculated pressures on the two prediction planes at the two time steps are shown in Figures 6  and 7. e calculated time signals at the two points on each plane are compared with the theoretical ones, shown in Figure 8.
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 P 1 are slightly better than those at the point P 2 . ese 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 T 1 and T 2 are defined by  T 1 � 〈p r x, y, z r , t p s x, y, z r , t 〉 �������������������������� 〈p 2 r x, y, z r , t 〉〈p 2 s x, y, z r , t 〉 , where 〈·〉 denotes the time averaged value, p s (x, y, z r , t) is the calculated pressure by the proposed method, and p r (x, y, z r , t) is the theoretical pressure. T 1 is a correlation coefficient which is sensitive to the similarity between the shapes of the signals and thus indicates their phase difference. e best value for T 1 would be 1. T 2 is the relative difference between the root mean square values of the two signals for characterizing the similarity of amplitudes. e target value to reach for T 2 is 0. e spatial maps of the indicators are shown in Figure 9. From the results, it can be seen that the values of indicator T 1 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 T 2 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. Furthermore, the indicators are also calculated at other Mach numbers on plane z r �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 T 1 at different Mach numbers are close to 1, more than 0.9 for all points. Meanwhile, the values of indicator T 2 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. e results indicate that the prediction errors are not amplified by the increasing flow effects as the velocity of the moving medium increases. is 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.

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][23][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
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.