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

Authors: Bichun Dong; Runmei Zhang; Chuanyang Yu; Huan Li
Journal: Shock and Vibration
 (2021)

---

## 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.


---

## Body

## 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 Section2. 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](1)1c2∂∂t+V⋅∇2px,y,z,t=∇2px,y,z,t,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=∂2/∂x2+∂2/∂y2+∂2/∂z2; ∇ is defined as ∇=∂/∂xi+∂/∂yj+∂/∂zk, 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 positivex direction; then, (1) can be written as(2)1c2∂2∂t2+2V∂2∂t∂x+V2∂2∂x2px,y,z,t=∇2px,y,z,t,where V is the magnitude of V.Taking the 2D spatial Fourier transform with respect tox and y, (2) can be expressed as(3)∂2∂t2+2iVkx∂∂t−V2kx2Pz,t=−c2kx2+ky2Pz,t+c2∂2Pz,t∂z2,where kx and ky represent the wavenumbers in x and y directions, respectively; i represents the imaginary unit; P represents the complex sound pressure in the wavenumber-time domain (here note that kx and ky 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. 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 thez direction can be represented as(4a)∂Pz,t∂t=Pz,t+Δt−Pz,t−Δt2Δt+OΔt2,(4b)∂2Pz,t∂t2=Pz,t+Δt−2Pz,t+Pz,t−ΔtΔt2+OΔt2,(4c)∂2Pz,t∂z2=Pz+Δz,t−2Pz,t+Pz−Δz,tΔz2+OΔz2,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.The substitution of (4) into (3) yields(5)1+iVkxΔtPz,t+Δt=r2Pz+Δz,t+Pz−Δz,t+21−r2−2b2Pz,t−1−iVkxΔtPz,t−Δt,where r is the grid ratio defined as r=cΔt/Δz, b is defined as b=cΔtKM/2 with KM=1−M2kx2+ky2, 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(6)1+iVkxΔtPm,n+1=r2Pm+1,n+Pm−1,n+21−r2−2b2Pm,n−1−iVkxΔtPm,n−1,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 wavenumber-time 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. The grid for calculation can be visually illustrated in Figure 1. To start the recursive calculation, it is given that P0,n represents the sound pressure on the initial plane to be propagated, and Pm,0 and Pm,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 P0,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 kx,ky.Figure 1
The grid of the recursive process in the wavenumber-time domain.The geometry of the proposed method is shown in Figure2, 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.Figure 2
The geometry of the FDWTD method.


### 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 Pkz,n of pressure Pm,n can be obtained from the spatial discrete Fourier transform:(7)Pkz,n=∑mPm,neikzmΔz,where kz is the wavenumber in z direction.It can be obtained from (7) that the wavenumber spectrum P±lkz,n of pressure Pm±l,n shifted by the spatial steps ±l is able to be represented as(8)P±lkz,n=Pkz,ne±ikzlΔz.The substitution of (7) and (8) into (6) yields(9)1+iVkxΔtPkz,n+1=2aPkz,n−1−iVkxΔtPkz,n−1,where a=1−r2−2b2+r2coskzΔz.Taking the Z-transform with respect to time in (9), it can be represented as(10)1+iVkxΔtZPkz,Z=2aPkz,Z−1−iVkxΔtZ−1Pkz,Z,in which Z denotes the variable after Z-transform, and thus the characteristic or amplification equation can be obtained as(11)1+iVkxΔtZ2−2aZ+1−iVkxΔt=0.The stability condition of (6) is that, for all values of kz 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(12)Z1,2=a±a2−1+V2kx2Δt21+iVkxΔt≤1.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 a2 and 1+V2kx2Δt2. It can be analyzed briefly as follows.Case 1.
a2>1+V2kx2Δt2
When the sign in inequality (12) takes “+,” the expression can be obtained as(13)a+a2−1+V2kx2Δt22≤1+V2kx2Δt2.
It can be obviously found that inequality (13) could not be satisfied for all values of kz.
When the sign in inequality (12) takes “–,” it yields(14)a2≥1+V2kx2Δt2,which cannot be satisfied for each wavenumber kz in the case of taking equality in inequality (14).Case 2.
a2≤1+V2kx2Δt2
Under this condition, whether the sign in inequality (12) takes “+” or “–,” it yields the same expression as follows:(15)1+V2kx2Δt2≤1+V2kx2Δt2.
This inequality can be always satisfied unconditionally.
Based on the analysis above, it can be concluded thatZ1,2≤1 can be unconditionally satisfied when a2≤1+V2kx2Δt2, while Z1,2≤1 cannot be always satisfied when a2>1+V2kx2Δt2. Thus, the expression of stability condition constrained by inequality (12) is(16)a2≤1+V2kx2Δt2.
Substitutinga=1−r2−2b2+r2coskzΔz into inequality (16) and applying the trigonometric formula, it yields(17)1−2b2−2r2sin2kzΔz22≤1+V2kx2Δt2.
This condition should be satisfied for all values ofkz, and thus it can be rewritten as(18)1−2b2−2r22≤1+V2kx2Δt2.
Substitutingb=cΔtKM/2 into inequality (18), the stability condition eventually can be represented as(19)r≤1+1+V2kx2Δt22−πFtfs2,where Ft is defined as Ft=cKM/2π and fs 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 ofPkz,n can be obtained as(20)Pkz,ω=∑nPkz,ne−iωnΔt,where ω represents the angular frequency.It can be obtained from (20) that the frequency spectrum P±gkz,ω of Pkz,n±g shifted by the time steps ±g is able to be represented as(21)P±gkz,ω=Pkz,ωe∓iωgΔt.Substituting (20) and (21) into (9), it yields(22)a=12e−iωΔt+eiωΔt+iVkxΔte−iωΔt−iVkxΔteiωΔt=cosωΔt+VkxΔtsinωΔt.Here it is worthy to note that, from the definition following (9), a is an even function of the wavenumber kx, and thus (22) should be actually represented as(23)a=cosωΔt+VkxΔtsinωΔt=1+V2kx2Δt2sinωΔt+φ,where the sign “⋅” represents the absolute value and φ is defined as φ=arctan1/VkxΔt.It has shown thata2 should satisfy the condition a2≤1+V2kx2Δt2 in the stability analysis. Now it is important to choose a proper specific value for a2, 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 thata2=1+V2kx2Δt2 and according to inequality (18), the relationship between b and r can be obtained as(24)b2+r2=1+1+V2kx2Δt22.Combining (23) and (24) and the definition of a, the dispersion relationship between the angular frequency ω and the wavenumber kz can be obtained as(25)1+V2kx2Δt2sinωΔt+φ=1−2b2−21+1+V2kx2Δt22−b2sin2kzΔz2.Substitutingb=πFt/fs and Δz=cΔt/r into (25), the wavenumber kz in a uniformly moving medium normalized by k0 can be yielded as(26)kzk0=1πf/fs1+1+V2kx2Δt22−πFtfs2×arcsin1−2πFt/fs2−1+V2kx2Δt2sin2πfΔt+φ1+1+V2kx2Δt2−2πFt/fs2,in which k0 is the acoustic wavenumber in the stationary medium defined as k0=ω/c=2πf/c.It can be found that there exists a transition frequencyft at which kz/k0 changes from the imaginary value to the real value, and the frequency ft normalized by the sampling frequency fs can be represented as(27)ftfs=12πarcsin11+V2kx2Δt21−2πFtfs2−φ.In fact, the theoretical value of wavenumberk˜z in the continuous case normalized by k0 in the propagation direction can be represented as(28)k˜zk0=1−1/πfsf/fsVkx+Ft/fs2f/fs2,where k˜z=k02−2k0Mkx−1−M2kx2−ky2 [5]. From (28), the normalized transition frequency f˜t/fs can be obtained as(29)f˜tfs=Vkx2πfs+Vkx2πfs2+Ftfs2.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 kx,ky, and the sampling frequency. The deviation indicates that, in the normalized frequency band between f˜t/fs and ft/fs, 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 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 ft/fs equals f˜t/fs, and thus the corrected value of Ft/fs satisfying the assumption can be obtained as(30)F˜tfs=1π121−1+V2kx2Δt2sinVkxfs+Vkxfs2+2πFtfs2+φ.Figure 3
Imaginary and real parts ofkz/k0 in the continuous and discretized cases at two different wavenumber points: kx,ky= (20, 20) [(a) and (b)] and kx,ky= (40, 40) [(c) and (d)].
(a)(b)(c)(d)The comparison at the wavenumber pointkx,ky= (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.Figure 4
Imaginary part (a) and real part (b) ofkz/k0 in the continuous case and those in the discretized case after correction at the wavenumber point kx,ky= (40, 40).
(a)(b)So far, the case wherea2=1+V2kx2Δt2 has been considered. If a2<1+V2kx2Δt2, the relationship between b and r becomes b2+r2<1/21+1+V2kx2Δt2, and a corresponding analysis can be conducted similarly to the previous process. However, compared with the case in which a2=1+V2kx2Δt2, a degradation of the results for kz/k0 can be observed. For instance, when a2=1/4 and b2=3/4−r2 are chosen, the results at the wavenumber point kx,ky= (40, 40) before and after corrections are presented in Figure 5.Figure 5
Imaginary and real parts ofkz/k0 at the wavenumber point kx,ky= (40, 40): uncorrected values [(a) and (b)] and corrected values [(c) and (d)].
(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 ofkz/k0. 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 tor2=1/21+1+V2kx2Δt2−b2. Substituting this condition and replacing b with the corrected value b˜=πF˜t/fs, the recursive expression of the wavenumber-time domain method used from now on can be finally obtained as(31)Pm,n+1=11+iVkxΔt1+1+V2kx2Δt22−b˜2Pm+1,n+Pm−1,n+1−1+V2kx2Δt2−2b˜2Pm,n−1−iVkxΔtPm,n−1.


## 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](1)1c2∂∂t+V⋅∇2px,y,z,t=∇2px,y,z,t,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=∂2/∂x2+∂2/∂y2+∂2/∂z2; ∇ is defined as ∇=∂/∂xi+∂/∂yj+∂/∂zk, 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 positivex direction; then, (1) can be written as(2)1c2∂2∂t2+2V∂2∂t∂x+V2∂2∂x2px,y,z,t=∇2px,y,z,t,where V is the magnitude of V.Taking the 2D spatial Fourier transform with respect tox and y, (2) can be expressed as(3)∂2∂t2+2iVkx∂∂t−V2kx2Pz,t=−c2kx2+ky2Pz,t+c2∂2Pz,t∂z2,where kx and ky represent the wavenumbers in x and y directions, respectively; i represents the imaginary unit; P represents the complex sound pressure in the wavenumber-time domain (here note that kx and ky 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. 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 thez direction can be represented as(4a)∂Pz,t∂t=Pz,t+Δt−Pz,t−Δt2Δt+OΔt2,(4b)∂2Pz,t∂t2=Pz,t+Δt−2Pz,t+Pz,t−ΔtΔt2+OΔt2,(4c)∂2Pz,t∂z2=Pz+Δz,t−2Pz,t+Pz−Δz,tΔz2+OΔz2,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.The substitution of (4) into (3) yields(5)1+iVkxΔtPz,t+Δt=r2Pz+Δz,t+Pz−Δz,t+21−r2−2b2Pz,t−1−iVkxΔtPz,t−Δt,where r is the grid ratio defined as r=cΔt/Δz, b is defined as b=cΔtKM/2 with KM=1−M2kx2+ky2, 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(6)1+iVkxΔtPm,n+1=r2Pm+1,n+Pm−1,n+21−r2−2b2Pm,n−1−iVkxΔtPm,n−1,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 wavenumber-time 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. The grid for calculation can be visually illustrated in Figure 1. To start the recursive calculation, it is given that P0,n represents the sound pressure on the initial plane to be propagated, and Pm,0 and Pm,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 P0,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 kx,ky.Figure 1
The grid of the recursive process in the wavenumber-time domain.The geometry of the proposed method is shown in Figure2, 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.Figure 2
The geometry of the FDWTD method.


## 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 Pkz,n of pressure Pm,n can be obtained from the spatial discrete Fourier transform:(7)Pkz,n=∑mPm,neikzmΔz,where kz is the wavenumber in z direction.It can be obtained from (7) that the wavenumber spectrum P±lkz,n of pressure Pm±l,n shifted by the spatial steps ±l is able to be represented as(8)P±lkz,n=Pkz,ne±ikzlΔz.The substitution of (7) and (8) into (6) yields(9)1+iVkxΔtPkz,n+1=2aPkz,n−1−iVkxΔtPkz,n−1,where a=1−r2−2b2+r2coskzΔz.Taking the Z-transform with respect to time in (9), it can be represented as(10)1+iVkxΔtZPkz,Z=2aPkz,Z−1−iVkxΔtZ−1Pkz,Z,in which Z denotes the variable after Z-transform, and thus the characteristic or amplification equation can be obtained as(11)1+iVkxΔtZ2−2aZ+1−iVkxΔt=0.The stability condition of (6) is that, for all values of kz 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(12)Z1,2=a±a2−1+V2kx2Δt21+iVkxΔt≤1.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 a2 and 1+V2kx2Δt2. It can be analyzed briefly as follows.Case 1.
a2>1+V2kx2Δt2
When the sign in inequality (12) takes “+,” the expression can be obtained as(13)a+a2−1+V2kx2Δt22≤1+V2kx2Δt2.
It can be obviously found that inequality (13) could not be satisfied for all values of kz.
When the sign in inequality (12) takes “–,” it yields(14)a2≥1+V2kx2Δt2,which cannot be satisfied for each wavenumber kz in the case of taking equality in inequality (14).Case 2.
a2≤1+V2kx2Δt2
Under this condition, whether the sign in inequality (12) takes “+” or “–,” it yields the same expression as follows:(15)1+V2kx2Δt2≤1+V2kx2Δt2.
This inequality can be always satisfied unconditionally.
Based on the analysis above, it can be concluded thatZ1,2≤1 can be unconditionally satisfied when a2≤1+V2kx2Δt2, while Z1,2≤1 cannot be always satisfied when a2>1+V2kx2Δt2. Thus, the expression of stability condition constrained by inequality (12) is(16)a2≤1+V2kx2Δt2.
Substitutinga=1−r2−2b2+r2coskzΔz into inequality (16) and applying the trigonometric formula, it yields(17)1−2b2−2r2sin2kzΔz22≤1+V2kx2Δt2.
This condition should be satisfied for all values ofkz, and thus it can be rewritten as(18)1−2b2−2r22≤1+V2kx2Δt2.
Substitutingb=cΔtKM/2 into inequality (18), the stability condition eventually can be represented as(19)r≤1+1+V2kx2Δt22−πFtfs2,where Ft is defined as Ft=cKM/2π and fs 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 ofPkz,n can be obtained as(20)Pkz,ω=∑nPkz,ne−iωnΔt,where ω represents the angular frequency.It can be obtained from (20) that the frequency spectrum P±gkz,ω of Pkz,n±g shifted by the time steps ±g is able to be represented as(21)P±gkz,ω=Pkz,ωe∓iωgΔt.Substituting (20) and (21) into (9), it yields(22)a=12e−iωΔt+eiωΔt+iVkxΔte−iωΔt−iVkxΔteiωΔt=cosωΔt+VkxΔtsinωΔt.Here it is worthy to note that, from the definition following (9), a is an even function of the wavenumber kx, and thus (22) should be actually represented as(23)a=cosωΔt+VkxΔtsinωΔt=1+V2kx2Δt2sinωΔt+φ,where the sign “⋅” represents the absolute value and φ is defined as φ=arctan1/VkxΔt.It has shown thata2 should satisfy the condition a2≤1+V2kx2Δt2 in the stability analysis. Now it is important to choose a proper specific value for a2, 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 thata2=1+V2kx2Δt2 and according to inequality (18), the relationship between b and r can be obtained as(24)b2+r2=1+1+V2kx2Δt22.Combining (23) and (24) and the definition of a, the dispersion relationship between the angular frequency ω and the wavenumber kz can be obtained as(25)1+V2kx2Δt2sinωΔt+φ=1−2b2−21+1+V2kx2Δt22−b2sin2kzΔz2.Substitutingb=πFt/fs and Δz=cΔt/r into (25), the wavenumber kz in a uniformly moving medium normalized by k0 can be yielded as(26)kzk0=1πf/fs1+1+V2kx2Δt22−πFtfs2×arcsin1−2πFt/fs2−1+V2kx2Δt2sin2πfΔt+φ1+1+V2kx2Δt2−2πFt/fs2,in which k0 is the acoustic wavenumber in the stationary medium defined as k0=ω/c=2πf/c.It can be found that there exists a transition frequencyft at which kz/k0 changes from the imaginary value to the real value, and the frequency ft normalized by the sampling frequency fs can be represented as(27)ftfs=12πarcsin11+V2kx2Δt21−2πFtfs2−φ.In fact, the theoretical value of wavenumberk˜z in the continuous case normalized by k0 in the propagation direction can be represented as(28)k˜zk0=1−1/πfsf/fsVkx+Ft/fs2f/fs2,where k˜z=k02−2k0Mkx−1−M2kx2−ky2 [5]. From (28), the normalized transition frequency f˜t/fs can be obtained as(29)f˜tfs=Vkx2πfs+Vkx2πfs2+Ftfs2.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 kx,ky, and the sampling frequency. The deviation indicates that, in the normalized frequency band between f˜t/fs and ft/fs, 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 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 ft/fs equals f˜t/fs, and thus the corrected value of Ft/fs satisfying the assumption can be obtained as(30)F˜tfs=1π121−1+V2kx2Δt2sinVkxfs+Vkxfs2+2πFtfs2+φ.Figure 3
Imaginary and real parts ofkz/k0 in the continuous and discretized cases at two different wavenumber points: kx,ky= (20, 20) [(a) and (b)] and kx,ky= (40, 40) [(c) and (d)].
(a)(b)(c)(d)The comparison at the wavenumber pointkx,ky= (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.Figure 4
Imaginary part (a) and real part (b) ofkz/k0 in the continuous case and those in the discretized case after correction at the wavenumber point kx,ky= (40, 40).
(a)(b)So far, the case wherea2=1+V2kx2Δt2 has been considered. If a2<1+V2kx2Δt2, the relationship between b and r becomes b2+r2<1/21+1+V2kx2Δt2, and a corresponding analysis can be conducted similarly to the previous process. However, compared with the case in which a2=1+V2kx2Δt2, a degradation of the results for kz/k0 can be observed. For instance, when a2=1/4 and b2=3/4−r2 are chosen, the results at the wavenumber point kx,ky= (40, 40) before and after corrections are presented in Figure 5.Figure 5
Imaginary and real parts ofkz/k0 at the wavenumber point kx,ky= (40, 40): uncorrected values [(a) and (b)] and corrected values [(c) and (d)].
(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 ofkz/k0. 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 tor2=1/21+1+V2kx2Δt2−b2. Substituting this condition and replacing b with the corrected value b˜=πF˜t/fs, the recursive expression of the wavenumber-time domain method used from now on can be finally obtained as(31)Pm,n+1=11+iVkxΔt1+1+V2kx2Δt22−b˜2Pm+1,n+Pm−1,n+1−1+V2kx2Δt2−2b˜2Pm,n−1−iVkxΔtPm,n−1.


## 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 conditionzmax>c1−M2T, where zmax and T represent the extremities of the analyzed region and observation time, respectively. Furthermore, the numerical stability and dispersion conditions require that r2=1/21+1+V2kx2Δt2−b2, and thus the grid interval Δz can be determined as(32)Δz=cΔt1/21+1+V2kx2Δt2−b2.From (32), it is noted that Δz changes with the values of b and kx, which eventually depends on the wavenumber point kx,ky. Therefore, if the propagation distance is given as d, the number of steps Nk that is necessary to reach the propagation distance for a certain wavenumber point can be determined by Nk=d/Δz.


### 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 zh=0.05 m. Two prediction planes are located at zr=0.2 m and 0.5 m. The sound pressure at r=x,y,z can be calculated as [18](33)pr,t=14πR1RR1q′t−Vx−x0R12qt,where R and R1 are represented as(34)R=−Mx−x0+R11−M2,R1=x−x02+1−M2y−y02+z−z02,where x0,y0,z0 is the coordinate of sound source, q′t represents the derivative of qt with respect to time, and qt represents the strength of the source given by(35)qt=0.005e−5.26300t−R/c−0.52.The 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.


### 3.2. Results and Discussions

In the simulations, the velocity of the moving medium presented as the Mach number is first set to beM=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 Δ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 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.Figure 6
The theoretical and calculated pressures on the prediction planezr=0.2 m when M=0.6: the theoretical pressures at 150 Δt (a) and 250 Δt (c) and the calculated pressures at 150 Δt (b) and 250 Δt (d).
(a)(b)(c)(d)Figure 7
The theoretical and calculated pressures on the prediction planezr=0.5 m when M=0.6: the theoretical pressures at 150 Δt (a) and 250 Δt (c) and the calculated pressures at 150 Δt (b) and 250 Δt (d).
(a)(b)(c)(d)Figure 8
The calculated time signals at the two pointsP1 and P2 on different prediction planes when M= 0.6: P1 (a) and P2 (b) on plane zr=0.2 m and P1 (c) and P2 (d) on plane zr=0.5 m.
(a)(b)(c)(d)It can be seen from Figures6 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 indicatorsT1 and T2 are defined by(36)T1=prx,y,zr,tpsx,y,zr,tpr2x,y,zr,tps2x,y,zr,t,T2=pr2x,y,zr,t−ps2x,y,zr,tpr2x,y,zr,t,where ⋅ denotes the time averaged value, psx,y,zr,t is the calculated pressure by the proposed method, and prx,y,zr,t 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.Figure 9
The spatial maps of indicatorsT1 and T2 to assess the phase and the amplitude similarities between the calculated time signals and the theoretical time signals when M=0.6 on plane zr=0.2 m [(a) and (b)] and on plane zr=0.5 m [(c) and (d)].
(a)(b)(c)(d)Furthermore, the indicators are also calculated at other Mach numbers on planezr=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.Figure 10
The temporal indicatorsT1 (a) and T2 (b) to assess the phase and the amplitude similarities between the calculated time signals and the theoretical time signals on plane zr=0.2 m at different Mach numbers.
(a)(b)


## 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 zh=0.05 m. Two prediction planes are located at zr=0.2 m and 0.5 m. The sound pressure at r=x,y,z can be calculated as [18](33)pr,t=14πR1RR1q′t−Vx−x0R12qt,where R and R1 are represented as(34)R=−Mx−x0+R11−M2,R1=x−x02+1−M2y−y02+z−z02,where x0,y0,z0 is the coordinate of sound source, q′t represents the derivative of qt with respect to time, and qt represents the strength of the source given by(35)qt=0.005e−5.26300t−R/c−0.52.The 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.


## 3.2. Results and Discussions

In the simulations, the velocity of the moving medium presented as the Mach number is first set to beM=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 Δ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 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.Figure 6
The theoretical and calculated pressures on the prediction planezr=0.2 m when M=0.6: the theoretical pressures at 150 Δt (a) and 250 Δt (c) and the calculated pressures at 150 Δt (b) and 250 Δt (d).
(a)(b)(c)(d)Figure 7
The theoretical and calculated pressures on the prediction planezr=0.5 m when M=0.6: the theoretical pressures at 150 Δt (a) and 250 Δt (c) and the calculated pressures at 150 Δt (b) and 250 Δt (d).
(a)(b)(c)(d)Figure 8
The calculated time signals at the two pointsP1 and P2 on different prediction planes when M= 0.6: P1 (a) and P2 (b) on plane zr=0.2 m and P1 (c) and P2 (d) on plane zr=0.5 m.
(a)(b)(c)(d)It can be seen from Figures6 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 indicatorsT1 and T2 are defined by(36)T1=prx,y,zr,tpsx,y,zr,tpr2x,y,zr,tps2x,y,zr,t,T2=pr2x,y,zr,t−ps2x,y,zr,tpr2x,y,zr,t,where ⋅ denotes the time averaged value, psx,y,zr,t is the calculated pressure by the proposed method, and prx,y,zr,t 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.Figure 9
The spatial maps of indicatorsT1 and T2 to assess the phase and the amplitude similarities between the calculated time signals and the theoretical time signals when M=0.6 on plane zr=0.2 m [(a) and (b)] and on plane zr=0.5 m [(c) and (d)].
(a)(b)(c)(d)Furthermore, the indicators are also calculated at other Mach numbers on planezr=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.Figure 10
The temporal indicatorsT1 (a) and T2 (b) to assess the phase and the amplitude similarities between the calculated time signals and the theoretical time signals on plane zr=0.2 m at different Mach numbers.
(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.



---

*Source: 1637338-2021-11-11.xml*