Time Domain Simulation of Sound Waves Using Smoothed Particle Hydrodynamics Algorithm with Artificial Viscosity

Smoothed particle hydrodynamics (SPH), as a Lagrangian, meshfree method, is supposed to be useful in solving acoustic problems, such as combustion noise, bubble acoustics, etc., and has been gradually used in sound wave computation. However, unphysical oscillations in the sound wave simulation cannot be ignored. In this paper, an artificial viscosity term is added into the standard SPH algorithm used for solving linearized acoustic wave equations. SPH algorithms with or without artificial viscosity are both built to compute sound propagation and interference in the time domain. Then, the effects of the smoothing kernel function, particle spacing and Courant number on the SPH algorithms of sound waves are discussed. After comparing SPH simulation results with theoretical solutions, it is shown that the result of the SPH algorithm with the artificial viscosity term added attains good agreement with the theoretical solution by effectively reducing unphysical oscillations. In addition, suitable computational parameters of SPH algorithms are proposed through analyzing the sound pressure errors for simulating sound waves.


Introduction
In 1977, Lucy [1] and Gingold and Monaghan [2] independently pioneered the smoothed particle hydrodynamics (SPH) method for modeling astrophysical phenomena in three-dimensional space.As a Lagrangian, meshfree method, the SPH method can handle fluid dynamic problems with complicated and time-variant domain topologies, large density ranges and object separation, as shown in recent reviews by Springel [3], Liu and Liu [4] and Monaghan [5].Like other meshless methods being used in computational acoustics [6,7], the SPH method has the potential to be used to solve complicated acoustic problems, such as combustion noise, bubble acoustics, sound propagation in multiphase flows, and so on [8].
Recently, Wolfe [9] and Hahn [10] directly solved the fluid dynamic equations by the SPH method to obtain pressure perturbations in the computational domain.Subsequently, Zhang et al. [8,11] proposed using the SPH method to solve linearized acoustic wave equations in the quiescent fluid.However, the appearance of unphysical oscillations leads to the instability of SPH simulation results.
In order to eliminate the unphysical oscillations phenomena, Von Neumann and Richtmyer [12] proposed the methodology of artificial viscosity for simulating one-dimensional shock waves.Landshoff [13] subsequently added a linear artificial viscosity term to the Von Neumann-Richtmyer artificial viscosity, which could further smooth the oscillations.On this basis, Monaghan and Gingold [14], Monaghan and Poinracic [15] and Monaghan [16] developed another type of artificial viscosity (Monaghan-type artificial viscosity) for simulating shocks.After this, other modifications to the Monaghan-type artificial viscosity have also been proposed [17], but the effects of the modified Monaghan-type artificial viscosity are still being investigated.In recent years, many research fields [18][19][20][21][22] have adopted various forms of artificial viscosity to capture shock or model a physical Navier-Stokes viscosity.So far, Monaghan-type artificial viscosity [14][15][16] is the most widely used artificial viscosity in the SPH simulations.Therefore, this paper focuses on introducing the Monaghan-type artificial viscosity [14][15][16] into the SPH simulation of sound propagation and interference.
The present paper is organized as follows.In Section 2, the SPH algorithm is built for solving the linearized acoustic wave equations, and a Monaghan-type artificial viscosity [14][15][16] term is added to the momentum equation.In Section 3, SPH algorithms with or without artificial viscosity are both used to simulate sound propagation and interference, and the effects of artificial viscosity are discussed.In Section 4, the numerical errors of sound pressure obtained by SPH algorithms with different kinds of smoothing kernel functions, different particle spacings and Courant numbers are calculated and analyzed.On this basis, suitable computational parameters of SPH algorithms for simulating sound waves are recommended.Finally, the conclusion of the paper is drawn in Section 5.

SPH Algorithm for Acoustic Simulation
Functions in the SPH method are written as a particle approximation form.The continuous SPH integral representation for f(r) can be written as: where f is a function of the vector r, Ω is the volume of the integral, W is the smoothing kernel function and h is the smoothing length.In the SPH convention, the kernel approximation operator is marked by the angle bracket < >.The particle approximation for the spatial derivative can be described as: 1 ( ) ( ) where ri and rj indicate the position of particles i and j separately, rij = ri − rj, rij is the distance between particle i and particle j, N is the number of particles in the computational domain, mj is the mass of particle j and ρj is the density of particle j.Wij = W(rij, h) is the smoothing kernel function.h is the smoothing length that defines the influence area of the smoothing function Wij.
In order to evaluate the effects of the kernel function, three different smoothing kernel functions are used to simulate sound propagation and interference, which are the Gaussian kernel function [2], the cubic spline kernel function [23] and the quintic spline kernel function [24].
The Gaussian kernel function can be written as: where q = rij/h, D α is 1/(π 1/2 h), 1/(πh 2 ) and 1/(π 3/2 h 3 ) in one-, two-and three-dimensional space, respectively.As discussed in [25], the Gaussian kernel function is sufficiently smooth even for high orders of derivatives, and it is very stable and accurate, especially for disordered particles.However, it is not really compact, as it never goes to zero theoretically, unless q approaches infinity.
As described in Crespo's thesis [25,26], so far, the cubic spline kernel has been the most widely-used smoothing function in the SPH literature, since it resembles a Gaussian function while having a narrower compact support (it is equal to zero for q > 2).However, the second derivative of the cubic spline is a piecewise linear function, so the stability properties are inferior compared to other smoothing kernel functions.
The quintic spline kernel function is: where q = rij/h, D α is 1/(120 h), 7/(478 πh 2 ) and 3/(359 πh 3 ) in one-, two-and three-dimensional space, respectively.According to the discussion on the quintic spline kernel function by Cao et al. [27], the quintic spline kernel function is a high order spline function and is not sensitive to the disorder particle distribution.The compact support domain of the quintic spline kernel function is narrower than the Gaussian kernel function.
In this paper, the medium of sound propagation and interference is ideal gas; the process of sound propagation and interference is adiabatic; and the particle velocity u, sound pressure δp and the density change of particle δρ satisfy the following equations: where P0 denotes the quiescent pressure of a medium, ρ0 is the quiescent density of a medium and c0 represents the initial sound speed.
The particle approximation of the continuity equation can be written as: ( ) The standard SPH formulation of the momentum equation can be written as: where the L i p δ and L j p δ are the sound pressure of particles i and j, respectively.
The particle approximation of state equation for ideal gas can be described as: If adding the artificial viscosity term, the momentum equation is obtained as: where Пij is the Monaghan-type artificial viscosity [14] and is given by: where αП, βП are adjustable non-dimensional constant.
( ) . ci and cj separately indicate the sound speed of particle i and j, and both ci and cj are equal to c0.In Equation ( 12), the expression of ij φ can be written as: where the factor φ = 0.1hij, hij = (hi + hj)/2.hi and hj indicate the smoothing length of particles i and j, respectively, and both of them are equal to h.The smoothing length is constant in this work.
Algorithm 1 shows the detail steps of acoustic wave simulation.In Algorithm 1, the second order leap-frog integration [28] is used to compute L i δρ , L i u , since there is low memory storage required and high efficiency in the computation.The expressions can be written as: where t Δ is the time step.In practice, time step t Δ should be larger than that estimated using the Courant-Friedrichs-Levy (CFL) condition [25].In this work, this condition is represented by:

Sound Propagation and Interference Model
We first build a one-dimensional sound propagation model to confirm agreement between SPH simulation results and theoretical solutions.In this model, the sound propagates from x < 0 to x ≥ 0 (x denotes the direction of propagation), and the SPH computational region is from −10 m to 90 m.The sound source is a plane wave, and the region of sound source is from −10 m to 0 m, while the region of sound propagation is from 0 m to 90 m.The sound wave is at [−10, 0] when t = 0 s.When t = 0.25 s, the distance of sound propagation is 85 m, that is the sound wave is located at [−10, 85].
The sound pressure of the acoustic wave in x < 0, namely the source of the acoustic wave, can be written as: where t denotes time, x is the geometric position in the propagation direction, PA denotes the amplitude of the sound wave, ω is the angular frequency of the sound wave and k = ω/c0 is the wave number.The initial values of the amplitude and angular frequency are set to be PA = 50 Pa and ω = 50 rad/s, respectively.In addition, the sound speed c0 is set to be 340 m/s, and the density of the propagation medium is 1.0 kg/m 3 .Similarly, a sound interference model is built, as well.In the interference model, the region of sound interference goes from −10 m to 80 m.In the computational region, a sound wave with sound pressure of 30 Pa and angular frequency of 100 rad/s transmits along the x axis; at the same time, another sound wave with sound pressure of 10 Pa and angular frequency of 100 rad/s comes from the opposite direction.After 0.15 s, the sound interference is generated by two sound waves.
The SPH algorithms with and without added artificial viscosity are used to simulate the sound propagation and interference.These simulation results are then compared to theoretical solutions.

Numerical Simulations and Model Verification
Table 1 lists the computational parameters that are used in the fundamental SPH algorithms of sound propagation and interference.The artificial viscosity term (Equation ( 12)) is added to the fundamental SPH algorithms to simulate sound propagation and interference.
With the purpose of validating the accuracy of our algorithm, the SPH results with and without artificial viscosity are compared to the theoretical solution in Figure 1.The SPH results with artificial viscosity were calculated on the condition that the αП and βП are set to be 0.06 and 0.50, respectively.As shown in Figure 1, the theoretical solution is plotted using a solid line, and the SPH results are represented by scattered points.The scattered points are plotted at intervals of 20 for clearly identifying the SPH results.It can be observed that the two SPH simulation results both possess similar smoothing trends compared to the theoretical solution.Several peaks and valleys appear in the graph between 0 m and 80 m, and the sound pressure amplitudes of all three computational cases agree well with each other.Both the SPH simulation results (with or without artificial viscosity) are in good agreement with the theoretical solution.
In order to provide a better comparison, three magnified portions of Figure 1 are provided in Figure 2a-c, separately.Figure 2a,b describes the first and second peaks, respectively, showing the portions of the peaks bounded by the black frames in Figure 1.Since the two valleys in Figure 1 are nearly identical, only one of them is shown in detail in Figure 2c.
The SPH simulation result in the region 26 ≤ x ≤ 38 m are plotted at intervals of four in Figure 2a.From the figure, it can be seen that the sound pressure of theoretical solutions firstly increases and then decreases with increasing x in a smoothly changing wave.However, the standard SPH algorithm gives a different result, as can be seen in the black frame in the middle of the figure.It is clear that the sound pressure obtained with the standard SPH algorithm fluctuates dramatically around the theoretical value, where two unphysical peaks and one unphysical valley can be seen.In comparison, the SPH algorithm with artificial viscosity obtains a much better result, that is the simulation results optimized by the artificial viscosity are closer to the theoretical solution.Furthermore, all of the points are stable with this result, and no obvious unphysical peak values are obtained.Similarly, the simulation results in Figure 2b altered with artificial viscosity are better than those obtained using the standard SPH algorithm, as can be seen from the values in the two black frames.Therefore, both Figure 2a and Figure 2b show that the SPH algorithm with artificial viscosity can reduce unphysical oscillations in sound propagation simulation.In Figure 2c, the scattered points of the SPH results in the region between 48 m and 58 m are plotted at intervals of four.From the figure, it can also be seen that the SPH algorithm with artificial viscosity can obtain a smoother result by effectively reducing unphysical oscillations.
Figure 3 shows the changes of sound pressure during sound interference.The sound interference results are calculated on the condition that the αП and βП are set to 0.06 and 0.50, respectively.Figure 3a gives the sources of acoustic waves in a particle form.In Figure 3b-d, the two acoustic waves propagate along two opposite directions, and there is no interference.With time increasing, the two acoustic waves meet and generate interference, as shown in Figure 3e,f.Figure 4 presents a comparison among the SPH results of sound interference (with and without artificial viscosity) and the theoretical solution at time t = 0.15 s.In Figure 4, at time t = 0.15 s, the SPH results of sound interference (with and without artificial viscosity) are compared to the theoretical solution.It can be observed that both the SPH results of sound interference (with and without artificial viscosity) are in good agreement with the theoretical solution.

Results of Numerical Errors and Discussion
In order to evaluate the SPH simulation results, the SPH numerical error of the sound pressure ɛpre is calculated.The error is defined by: where ɛpre is a non-dimensional parameter and n denotes the number of particles that are used to compute numerical error.L j p δ is the simulation sound pressure of the particle j, and ( ) j p r δ is the theoretical sound pressure in the position of particle j.P denotes the pressure that is used to non-dimensionalize the numerical error.In this paper, P is defined as the amplitude of an acoustic wave in the sound propagation and is defined as the difference of the sound pressure amplitudes of two acoustic waves in the sound interference, that is P = 50 Pa in the sound propagation and P = 20 Pa in the sound interference.The region that is used to compute the numerical error is from 0 m to 90 m in the sound propagation and from 0 m to 70 m in the sound interference.
According to the computational parameters given in Table 1, we simulate the sound propagation and interference using the SPH algorithms with αП changing from 0.00 to 0.08; meanwhile, βП is set as 0.50.Then, we calculate the SPH numerical errors of the sound pressure according to Equation (18), as shown in Table 2.
From Table 2, it can be seen that the SPH algorithms with artificial viscosity can improve the accuracy of the simulation.Compared to the sound interference simulation, when the parameters of artificial viscosity term are the same, the SPH numerical errors of sound propagation are smaller.For sound propagation, the SPH numerical errors of sound pressure reach the minimum when αП = 0.06, βП = 0.5.For sound interference, the SPH numerical errors of sound pressure reach the minimum when αП = 0.03, βП = 0.5.Since sound propagation and interference are similar, we mainly discuss the effects of artificial viscosity on different sound interference models in the present work.In the discussion, different computational parameters, namely the smoothing kernel function, particle spacing and the Courant number, are considered.
Different sound interference cases with Gaussian, cubic spline and quintic spline kernel functions are simulated.The artificial viscosity terms with different values of αП and βП are added, and the numerical errors of sound pressure are plotted in Figure 5.As shown in Figure 5, for the same particle spacing Δx and Courant number, the smoothing kernel functions have obvious effects on the SPH simulation results.The SPH simulation using the cubic spline kernel function has more accuracy results than that using the Gaussian kernel function or the quintic spline kernel function.The numerical errors of SPH simulation using the quintic spline kernel function are slightly less than those using the Gaussian kernel function, but there is not much difference between these two types.
From Figure 5a, it can be found that, for different smoothing kernel functions, adding an artificial viscosity term can reduce the numerical errors of sound pressure.The numerical errors of the SPH results with three kinds of smoothing kernel functions share a similar changing tendency with the increasing of αП.They firstly decrease and then increase and reach their minimums when αП locates at [0.03, 0.06].However, as shown in Figure 5b, parameter βП has no obvious effects on the simulation results.
The sound interference models with particle spacing ∆x separately set to 0.04 m, 0.06 m and 0.08 m are also simulated.Figure 6 shows the numerical errors of sound pressure obtained by different sound interference models.As can be seen from Figure 6, for the same smoothing kernel function and Courant number, ignoring the effect of artificial viscosity terms, the largest numerical error of sound pressure occurs when particle spacing is 0.08 m, and the least error of sound pressure happens at a particle spacing of 0.04 m.That is, the numerical errors are heavily influenced by the particle spacing Δx, and the numerical errors increase with the particle spacing increasing.
Similarly, it can also be found that the numerical errors of the SPH results with different particle spacings all have the same changing tendency with the increasing of αП.Parameter αП has obvious effects on the simulation results, while parameter βП does not.
In addition, different sound interference models with Courant numbers separately set to 0.02, 0.04, 0.06 and 0.08 are simulated.Figure 7 shows the numerical errors of sound pressure for sound interference models with different Courant numbers.Figure 7 indicates that a better simulation result can be found when the Courant number is 0.08.For different Courant numbers, the numerical errors of sound pressure still share the same changing tendency with parameter αП increasing.Similarly, parameter βП has no obvious effect on the simulation results.
In summary, the SPH algorithms using the cubic spline have optimum simulation results when the particle spacing ∆x = 0.04 m and the Courant number is 0.08.The corresponding numerical errors could reach their minimum value when αП belongs to the region [0.03, 0.06].

Conclusions
In this work, an artificial viscosity term is added to the standard SPH algorithm in order to improve its accuracy in modelling sound propagation and interference.SPH forms of the continuity equation, the momentum equation and the equation of state are established firstly.Following this, a one-dimensional sound propagation model and a one-dimensional sound interference model are computed to evaluate the effects of the artificial viscosity term.A comparison between SPH results with and without artificial viscosity demonstrates that a more stable and smoother simulation result can be obtained by adding an artificial viscosity term to the standard SPH algorithm, which is due to effectively reducing the unphysical oscillations.Furthermore, analysis of numerical errors of sound pressure shows that the SPH simulation using the cubic spline kernel function can obtain more accurate results with the particle spacing being 0.04 m and the Courant number being 0.08.Parameter αП has more obvious effects than βП on the simulation results, and the SPH algorithms of sound waves can obtain optimal results when αП is in the region between 0.03 and 0.06, considering sound pressure error.

u
where L i δρ represents the density change of particle i, L i ρ and L j ρ separately denote the density of particles i and j, the superscript L stands for the Lagrangian variable associated with a fluid particle, t is the time, L j m is the mass of particle j, are the velocity of particles i and j, respectively.

Algorithm 1 .δρ
The SPH algorithmic description for acoustic wave simulation.begin Initialization phase Generate sound source according to function of sound pressure ( ) p = f t,x while (time step < total time step) for all particle i do Calculate the spatial derivative of smoothing kernel i ij W ∇ Compute the derivative of the density change with respect to time ( ) Compute the derivative of the velocity with respect to time L i D Dt u Introduce the artificial viscosity term, and recompute L i D Dt u Compute L i δρ , L i u using the second order leap-frog integration end for Recompute the sound source according to function of sound pressure ( ) p = f t,x end while Return the simulation result, which includes sound pressure, density change, and so on.end

Figure 1 .
Figure 1.Sound propagation pressure of acoustic wave at t = 0.25 s.

Figure 2 .
Figure 2. Comparison of SPH results (with and without artificial viscosity) with theoretical solutions.(a) Detailed view of the first peak; (b) detailed view of the second peak; (c) detailed view of the valley.

Figure 4 .
Figure 4. Sound interference pressure of acoustic waves at t = 0.15 s.

Figure 5 .
Figure 5. Numerical errors of sound pressure for sound interference models with different kinds of smoothing kernel functions.Particle spacing Δx is 0.04 m, and the Courant number is 0.04.(a) Effects of αП on the simulation results at βП = 0.5; (b) effects of βП on the simulation results at αП = 0.0.

Figure 6 .
Figure 6.Numerical errors of sound pressure for sound interference models with different particle spacings.The smoothing kernel function is cubic spline, and the Courant number is 0.04.(a) Effects of αП on the simulation results at βП = 0.5; (b) effects of βП on the simulation results at αП = 0.0.

Figure 7 .
Figure 7. Numerical errors of sound pressure for sound interference models with different Courant numbers.The particle spacing Δx is 0.04 m and the smoothing kernel function is cubic spline.(a) Effects of αП on the simulation results at βП = 0.5; (b) effects of βП on the simulation results at αП = 0.0.

Table 1 .
Parameters for the SPH algorithms of sound propagation and interference.

Table 2 .
The SPH numerical errors of sound pressure.