A Lagrangian approach for computational acoustics with particle-based method

Engineering Analysis with Boundary Elements 108 (2019) 459–471

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

A Lagrangian approach for computational acoustics with particle-based method Yongou Zhang a,b,c, Stefan G. Llewellyn Smith b, Tao Zhang a,∗, Tianyun Li a a

School of Naval Architecture and Ocean Engineering, Huazhong University of Science and Technology, Wuhan 430074, China Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla, CA 92093, USA c School of Transportation, Wuhan University of Technology, Wuhan 430063, China b

a r t i c l e

i n f o

Keywords: Lagrangian approach Computational acoustics Particle-based method Smoothed particle hydrodynamics Perturbation equation

a b s t r a c t Although Eulerian approaches are standard in computational acoustics, they are less effective for problems with moving boundaries and multiphase systems like bubble acoustics. In this paper, a Lagrangian approach to model sound propagation in moving fluid is presented and implemented numerically with particle-based methods. Fluid dynamic equations is divided into a set of hydrodynamic equations for the motion of fluid particles and perturbation equations for the acoustic quantities corresponding to each fluid particle. Then, the smoothed particle hydrodynamics (SPH) method and the corrective smoothed particle (CSP) method are introduced to solve the perturbation equations in Lagrangian form. A hybrid meshfree and finite-difference time-domain boundary treatment technique is utilized to represent acoustic boundaries. Finally, applications to modeling sound propagation in steady or unsteady fluids in motion are outlined, treating a number of different cases in one and two space dimensions. The Lagrangian approach shows good agreement with exact solutions. The comparison indicates that the CSP method exhibits accuracy convergence in cases with different background flow. Different acoustic boundary conditions are validated as being effective for benchmark problems in computational acoustics.

1. Introduction Computational acoustics is widely used in industrial applications to model devices such as fans, automobile manufactures, underwater vehicles, room acoustics, and so on. To deal with these kinds of problems, different numerical solvers based on Eulerian or Lagrangian approaches are developed. These numerical solvers are mesh-based or meshfree. By using Eulerian mesh-based and meshfree methods, various acoustic problems in engineering have been solved so far. Different Eulerian mesh-based numerical methods include the finite difference method (FDM) [1], the finite element method (FEM) [2], the boundary element method (BEM) [3], and other modified or coupled methods [4–9]. Eulerian meshfree methods has been used to solve the acoustic wave equation in quiescent media in the time domain or the Helmholtz equation in the frequency domain, as in the multiple-scale reproducing kernel particle method (RKPM) [10], the element-free Galerkin method (EFGM) [11], the method of fundamental solutions (MFS) [12], etc. [13–17]. However, as an Eulerian approach, they are probably not the most suitable methods for certain acoustic problems, such as transient acoustic problems with free surfaces, complex material interfaces, moving or deformable boundaries, and multiphase systems. Some specific examples

∗

are bubble acoustics [18], combustion noise [19], sound propagation in multiphase flows [20], etc. The present paper focuses on a pure Lagrangian meshfree particlebased computational acoustics (PCA) method. Compared to the Eulerian approach which solves for quantities at fixed locations in space, the Lagrangian approach uses individual particles which move through both space and time, and have their own physical properties, such as density, velocity, and pressure, to represent the dynamically evolving fluid flow. The flow is described by recording the time history of each fluid particle. The application of Lagrangian approaches in CFD has shown the advantages of such methods in solving problems with multiphase systems [21–23]. The PCA method uses the Lagrangian approach, so it is available to simulate the sound generation and propagation from the perspective of fluid particles. It has the following advantages: (i) numerical error accompanying the calculation of advection term is eliminated since the advection term is included in the Lagrangian derivative; (ii) the Lagrangian characteristic makes it easy to handle complex changes of computational domain shape and the problems with moving boundaries; (iii) the interface between different media can be naturally tracked by difference in particle density without using any additional algorithm such as the volume of fluid model; (iv) it is available to achieve

Corresponding author. E-mail addresses: [email protected] (Y. Zhang), [email protected] (S.G.L. Smith), [email protected] (T. Zhang), [email protected] (T. Li).

https://doi.org/10.1016/j.enganabound.2019.07.020 Received 20 February 2019; Received in revised form 4 July 2019; Accepted 29 July 2019 Available online 19 September 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

particle respectively; p0 , 𝜌0 , and v0 are the components in the absence of sound, while 𝛿p, 𝛿𝜌, 𝛿v, and 𝛿S are the perturbations of corresponding parameters due to acoustic waves. All these variables are supposed to be the physical properties of a fluid particle which can move through both space and time, rather than at a fixed position. The Lagrangian form of the fluid dynamics equations including gravity can be written as: the continuity equation

the derivative solution in a local support domain rather than the entire computational domain so that parallel computing can be easily realized. The introduction of PCA method into acoustic simulation enables these advantages to be fully utilized in the field of computational acoustics, and provides a Lagrangian approach for the numerical research of acoustic problems. However, in the field of computational acoustics, there have been almost no applications of Lagrangian approaches. It is possible to solve the Navier–Stokes (N–S) equations directly over all spatial and temporal scales in a Lagrangian approach [24,25], which can be viewed as a direct numerical simulation (DNS), but this Lagrangian DNS is as costly as an Eulerian DNS. Another way is to use a partial Lagrangian method [26,27] which computes flow field with fixed grids, and hence the advantages of the Lagrangian approach mentioned before cannot be used in the fluid flow computation. The purpose of our work is to devise a Lagrangian approach for computational acoustics, which computes the fluid particle motion and the sound waves based on the acoustic perturbation assumption [28,29]. In this way, this approach can consider sound generation and propagation with complex background flow or multiple media, and can be more efficient than the Lagrangian DNS method by ignoring scale disparities since the flow and acoustics are treated separately [30,31] as in hydrodynamic/acoustic splitting methods. In addition, we choose particlebased method, instead of traditional mesh-based methods to realize the computation numerically, because dependent variables correspond fluid particles, which can move with the fluid, instead of fixed points. In this paper, we have chosen the smoothed particle hydrodynamics (SPH) method and the corrective smoothed particle (CSP) method to implement and validate the Lagrangian approach. These two methods are representative of particle-based methods in general. The SPH method is a Lagrangian meshfree method first pioneered independently by Lucy [32] and Gingold and Monaghan [33] to solve astrophysical problems in 1977. It has been widely used in different fields as outlined in recent reviews [21–23,34]. Details about the standard SPH computation can be found in Liu and Liu’s book [35], and the CSP method can be found in [36,37]. To date, the Lagrangian method like the SPH method, the finite particle method [38,39], the moving particle semi-implicit method [40], etc. have shown its advantages in fluid dynamics simulations, and the Lagrangian particle tracking method has been widely used in modeling bubble motions [41], free surface [42], dust [43], and other problems. In the preliminary work, we have given a compact discussion concerns the numerical performance of the SPH method in solving acoustic wave equations in quiescent media [44–46]. The present paper is organized as follows. In Section 2, the Lagrangian form of perturbation equation are obtained. In Section 3, the SPH and CSP formulations for simulating acoustic waves in a moving flow is given based on the Lagrangian form of perturbation equation. In Section 4, applications of particle-based methods to different acoustic problems are tested separately. Several test cases are given in order to validate the numerical method in Section 5, while Section 6 summarizes the results of this work.

D𝜌 = −𝜌∇ ⋅ 𝒗, D𝑡

(1)

the momentum equation D𝒗 1 = − ∇ 𝑃 + 𝒈, D𝑡 𝜌

(2)

where the material or Lagrangian derivative is written as 𝐷 𝜕 = + 𝒗 ⋅ ∇. 𝐷𝑡 𝜕𝑡

(3)

Now we consider the sound waves in an inhomogeneous medium without the effects of viscosity. Sound propagation can be modeled as an adiabatic process, and viscosity, diffusion and thermal conductivity can be neglected. Let c be the sound velocity; then the equation of state becomes 𝐷𝜌 𝐷𝑃 = 𝑐2 , 𝐷𝑡 𝐷𝑡

(4)

Eqs. (1), (2), and (4) are a closed system for equations of fluid dynamics. 2.2. Equations for background fluid motion Now we consider the flow without sound waves. The usual hydrodynamic equations include gravity can be written as 𝐷 0 𝜌0 = − 𝜌0 ∇ ⋅ 𝒗0 , 𝐷𝑡 D0 𝒗0 1 = − ∇ 𝑝 0 + 𝒈, D𝑡 𝜌0

(5) (6)

where the material derivative is 𝐷0 𝜕 = + 𝒗0 ⋅ ∇. 𝐷𝑡 𝜕𝑡

(7)

For an adiabatic process neglecting viscosity, diffusion, and thermal conductivity, we have the equation of state 𝐷0 𝑝0 𝐷 𝜌 = 𝑐02 0 0 , 𝐷𝑡 𝐷𝑡

(8)

where c0 is the sound speed not influenced by the acoustic wave, with 𝑐 2 = 𝑐02 + 𝛿𝑐 2 . 2.3. Lagrangian form of perturbation equations

2. Lagrangian acoustic perturbation equations

Considering acoustic perturbations generated by sound, the continuity equation (Eq. (1)), the momentum equations (Eq. (2)), and the equation of state (Eq. (4)) lead to

The fundamental assumptions and equations are presented first, and then the Lagrangian form of perturbation equation for computational acoustics are obtained.

𝐷(𝜌0 + 𝛿𝜌) = −(𝜌0 + 𝛿𝜌)∇ ⋅ (𝒗0 + 𝛿𝒗), 𝐷𝑡 𝐷(𝒗0 + 𝛿𝒗) 1 =− ∇(𝑝0 + 𝛿𝑝) + 𝒈, 𝐷𝑡 (𝜌0 + 𝛿𝜌) 𝐷(𝑝0 + 𝛿𝑝) 𝐷(𝜌0 + 𝛿𝜌) = (𝑐02 + 𝛿𝑐 2 ) . 𝐷𝑡 𝐷𝑡

2.1. Fundamental assumptions and equations The basic idea of the hydrodynamic/acoustic splitting method is to split the background flow and acoustic perturbation based on differences in amplitude. Details of this approach can be found in [28,29]. Following this idea, consider a fluid particle which has physical properties: P = p0 + 𝛿p, |𝛿p| < < p0 ; 𝜌 = 𝜌0 + 𝛿𝜌, |𝛿𝜌| < <𝜌0 ; v = v0 + 𝛿v, |𝛿v| < < v0 ; P, 𝜌, and v are the pressure, density, and velocity of the fluid

(9) (10) (11)

Rearranging these equations using the hydrodynamic equations (Eqs. (5), (6), and (8)), we obtain the Lagrangian form of perturbation equations (PE) as 460

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

𝐷0 𝛿𝜌 = −𝜌0 ∇ ⋅ 𝛿𝒗 − 𝛿𝜌∇ ⋅ (𝒗0 + 𝛿𝒗) − (𝛿𝒗 ⋅ ∇)(𝜌0 + 𝛿𝜌), 𝐷𝑡

(12)

𝐷0 𝛿𝒗 𝛿𝜌 1 = − ∇𝛿𝑝 − (𝛿𝒗 ⋅ ∇)(𝒗0 + 𝛿𝒗) + ∇(𝑝0 + 𝛿𝑝), 𝐷𝑡 𝜌0 𝜌0 (𝜌0 + 𝛿𝜌)

(13)

𝐷0 𝛿𝑝 𝐷 𝛿𝜌 − 𝑐02 0 = −𝛿𝒗∇(𝑝0 + 𝛿𝑝) + 𝑐02 (𝛿𝒗 ⋅ ∇)(𝜌0 + 𝛿𝜌) 𝐷𝑡 𝐷𝑡 [ ] 𝐷 𝛿𝜌 2 𝐷 0 𝜌0 + 𝛿𝑐 + 0 + (𝛿𝒗 ⋅ ∇)(𝜌0 + 𝛿𝜌) , 𝐷𝑡 𝐷𝑡

where f is a function of the vector r, Ω is the volume of the integral, W is the smoothing kernel, and h is the smoothing length which represents the size of the support domain. The kernel approximation operator is denoted by the angle bracket <>. The particle approximation for the function f (r) at particle i can be written as < 𝑓 ( 𝒓𝑖 ) > =

(14)

𝐷0 𝛿𝑝 𝐷 𝛿𝜌 𝐷 𝜌 − 𝑐02 0 = 𝑐02 (𝛿𝒗 ⋅ ∇)𝜌0 − (𝛿𝒗 ⋅ ∇)𝑝0 + 𝛿𝑐 2 0 0 , 𝐷𝑡 𝐷𝑡 𝐷𝑡

< ∇ 𝑓 ( 𝒓𝑖 ) > =

(15)

(19)

𝑁 ∑ 𝑚𝑗 𝑗=1

𝜌𝑗

𝑓 (𝒓𝑗 )∇𝑖 𝑊𝑖𝑗 ,

𝒓𝑖 −𝒓𝑗 𝜕 𝑊𝑖𝑗

(20)

(17)

the divergence of a vector field. The cubic spline function is chosen as the smoothing kernel in the present paper.

𝑟𝑖𝑗

𝜕 𝑟𝑖𝑗

and rij = |ri − rj |. This can also be derived for

3.2. SPH equations for background flow Particle approximation for fluid dynamic equations (Eqs. (1), (2), and (4)) in the absence of forces are given in [35]. The particle approximation used for the hydrodynamic equations including gravity (Eqs. (5), (6), and (8)) can be obtained in a similar way. The results are

Since the SPH method is a famous and widely used particle-based method in computational fluid dynamics, the linearized Lagrangian form of perturbation equations are solved using SPH and a correction to it, CSP, in this section.

𝑁 ∑ 𝑚𝑗 𝐷 0 𝜌𝑖 = 𝜌𝑖 𝒗 ⋅ ∇𝑖 𝑊𝑖𝑗 , 𝐷𝑡 𝜌 𝑖𝑗 𝑗=1 𝑗

(21)

𝑁 ∑ 𝑝𝑖 + 𝑝𝑗 𝐷 0 𝒗𝑖 =− 𝑚𝑗 ∇𝑖 𝑊𝑖𝑗 + 𝒈, 𝐷𝑡 𝜌𝑖 𝜌𝑗 𝑗=1

(22)

𝐷0 𝑝𝑖 𝐷 𝜌 = 𝑐02 0 𝑖 , 𝐷𝑡 𝐷𝑡

(23)

where pi , 𝜌i , and vi are the pressure, density, and velocity of the fluid particle i in the absence of sound, corresponding to p0 , 𝜌0 , and v0 in Section 2, vij = vi − vj , and ∇i denotes the gradient of the kernel taken with respect to the coordinate of particle i. Other forms of (22) and (23) are

3.1. Basic concepts of SPH SPH is a meshfree Lagrangian particle method. Functions in the SPH method are represented by a particle approximation. Some basic concepts are presented in this section. A function f (r) can be represented as ∫Ω

𝑓 (𝒓𝑗 )𝑊𝑖𝑗 ,

where ∇𝑖 𝑊𝑖𝑗 =

3. SPH and CSP equations

𝑓 (𝒓′ )𝑊 (𝒓 − 𝒓′ , ℎ)𝑑 𝒓′ ,

𝜌𝑗

(16)

where 𝛿𝜌a represents a volume velocity generated by a source, and 𝛿f/𝜌0 is a force density [28]. This corresponds to a singular source in the fluid, which can be useful for modeling purposes. Here 𝛿c2 be the perturbation of the square of sound velocity caused by the sound wave. It can be seen that the background motion of the fluid particle can be obtained from Eqs. (5), (6), and (8), and the acoustic variables of each particle can be obtained from Eqs. (15)–(17).

< 𝑓 ( 𝒓) > =

𝑗=1

where ri and rj are the positions of particles i and j, N is the number of particles in the computational domain, mj is the mass of particle j, Wij = W(rij ,h) and rij is the distance between particle i and particle j as shown in Fig. 1. The spatial derivative ∇f(ri ) is similarly obtained as

where the material derivative used is Eq. (8). These equations are exact. Ignoring the second-order terms, leads to linear equations for sound waves. In order to take into account sound sources, two terms are added to the right hand side of equations, and the linearized Lagrangian form of perturbation equations (LPE) are 𝐷0 𝛿𝜌 = −𝜌0 ∇ ⋅ 𝛿𝒗 − 𝛿𝜌∇ ⋅ 𝒗0 − (𝛿𝒗 ⋅ ∇)𝜌0 + 𝛿𝜌𝑎, 𝐷𝑡 𝐷0 𝛿𝒗 𝛿𝒇 𝛿𝜌 1 = − ∇𝛿𝑝 − 𝛿𝒗∇ ⋅ 𝒗0 + ∇𝑝0 + , 𝐷𝑡 𝜌0 𝜌0 (𝜌0 )2

𝑁 ∑ 𝑚𝑗

𝜌𝑖 =

(18)

𝑁 ∑ 𝑗=1

𝑚𝑗 𝑊𝑖𝑗 ,

(24)

Fig. 1. Fluid particles, support domain and the kernel function.

461

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

[ ] 𝑁 ∑ 𝑝𝑗 𝐷0 𝒗𝑖 𝑝𝑖 =− 𝑚𝑗 + ∇𝑖 𝑊𝑖𝑗 + 𝒈. 𝐷𝑡 (𝜌𝑖 )2 (𝜌𝑗 )2 𝑗=1

+

(25)

1 𝛼 (𝒓 − 𝝃 𝛼 )(𝒓𝛽 − 𝝃 𝛽 )𝑓𝛼𝛽 (𝝃)𝑊 (𝒓−𝝃, ℎ)𝑑 𝒓 + .... ∫Ω 2!

(33)

Ignoring all derivative terms, f (𝝃) can be represented as 3.3. Perturbation equations in SPH form < 𝑓 (𝝃) >=

Applying the SPH particle approximation equation Eq. (20) to the linearized Lagrangian perturbation equations (Eqs. (15)–(17)) yields 𝑁 ∑

𝑚𝑗 𝑚𝑗 𝐷0 𝛿𝜌𝑖 = − 𝜌𝑖 𝛿𝒗𝑗 ⋅ ∇𝑖 𝑊𝑖𝑗 − 𝛿𝜌𝑖 𝒗 ⋅ ∇𝑖 𝑊𝑖𝑗 𝐷𝑡 𝜌 𝜌 𝑗 𝑗 𝑗=1 𝑗=1 𝑗 − 𝛿𝒗𝑖

𝑗=1

𝜌𝑗

𝜌𝑗 ∇𝑖 𝑊𝑖𝑗 + 𝛿𝜌𝑖 𝑎𝑖 ,

∑𝑁 𝑓 ( 𝒓𝑖 ) = (26) ( ) 𝑓𝛼 𝒓𝑖 =

𝑁 𝑁 ∑ 𝑚𝑗 𝐷0 𝛿𝒗𝑖 1 ∑ 𝑚𝑗 =− 𝛿𝑝 ∇ 𝑊 − 𝛿𝒗𝑖 ⋅ 𝒗 ⋅ ∇𝑖 𝑊𝑖𝑗 𝐷𝑡 𝜌𝑖 𝑗=1 𝜌𝑗 𝑗 𝑖 𝑖𝑗 𝜌 𝑗 𝑗=1 𝑗

+

𝑁 𝛿𝜌𝑖 ∑ 𝑚𝑗 𝛿𝒇 𝑖 𝑝𝑗 ∇𝑖 𝑊𝑖𝑗 + , 𝜌𝑖 (𝜌𝑖 )2 𝑗=1 𝜌𝑗

𝐷0 𝛿𝑝𝑖 𝐷 𝛿𝜌 − 𝑐02 0 𝑖 𝐷𝑡 𝐷𝑡 𝑁 𝑁 ∑ ∑ 𝑚𝑗 𝑚𝑗 𝐷 𝜌 = 𝑐02 𝛿𝒗𝑖 ⋅ 𝜌𝑗 ∇𝑖 𝑊𝑖𝑗 − 𝛿𝒗𝑖 ⋅ 𝑝𝑗 ∇𝑖 𝑊𝑖𝑗 + 𝛿𝑐 2 0 𝑖 , 𝜌 𝜌 𝐷𝑡 𝑗=1 𝑗 𝑗=1 𝑗

𝑚𝑗

𝑗=1

𝜌𝑗

𝜌𝑗 ∇𝑖 𝑊𝑖𝑗 + 𝛿𝜌𝑖 𝑎𝑖 ,

(27)

𝑁 𝛿𝜌𝑖 ∑ 𝑝𝑖 + 𝑝𝑗 𝛿𝒇 𝑖 𝑚 ∇𝑖 𝑊𝑖𝑗 + . 𝜌𝑖 𝑗=1 𝑗 𝜌𝑖 𝜌𝑗 𝜌𝑖

𝜕𝑓 (𝝃) , 𝜕 𝒓𝛼

𝑓𝛼𝛽 (𝝃) =

𝜕 2 𝑓 (𝝃) , 𝜕 𝒓𝛼 𝜕 𝒓𝛽

(30)

∫Ω

𝑓 (𝝃)𝑊 (𝒓−𝝃, ℎ)𝑑 𝒓 +

∫Ω

(35)

[ ( ) ( )] 𝑓 𝒓𝑗 − 𝑓 𝒓𝑖 ∇𝑖 𝑊𝑖𝑗 . ) ∑𝑁 𝑚𝑗 ( 𝛼 𝒓𝑗 − 𝒓𝛼𝑖 ∇𝑖 𝑊𝑖𝑗 𝑗=1 𝜌

∑𝑁

𝑚𝑗 𝑗=1 𝜌𝑗

(36)

if 𝑖 = f luid particles if 𝑖 = boundary particles if 𝑖 = virtual particles

where 𝛼 1 = 15/ c0 2 , 𝛼 2 = ln 2/16, c0 = 340 m/s, 𝜌0 is the initial density of medium, x represents the horizontal coordinates of the particle in the computation domain, while y represents the vertical coordinate of the particle. The simulation results of the Guassian pulse of sound propagation with rigid acoustic boundary are shown in Fig. 2. Initial particle spacing is set as 0.2 m and h = 0.22 m. Fig. 2(a), (b) gives the contour of sound pressure at t = 0.03 s and 0.09 s. It can be seen that the acoustic wave propagates to all around in the form of circle, and the center of circle is located at the origin of the coordinate. The speed of propagation is 340 m/s. As time goes on, the radius of the circle increases, while the amplitude of the sound waves decreases. Sound pressure and velocity perturbations of particles at y = 0 are compared with the theoretical solution, as shown in Fig. 3. In the figure, the red solid line represents the theoretical solution, the black scattered point represents the simulation results at y = 0 m.

(31)

(32)

𝑓 (𝒓)𝑊 (𝒓−𝝃, ℎ)𝑑 𝒓

=

,

To verify the feasibility and accuracy of two-dimensional sound propagation by using the CSPM codes, a Guassian pulse model given by ICASE/LaRC workshop is built, which are proposed for solving the benchmark problem in computational aeroacoustics (CAA). At the time t = 0, the distribution of the sound pressure of Guassian pulse model in the computation domain can be written as [ ] 𝛿𝑝(𝑡 = 0, 𝑥, 𝑦) = 𝛼1 𝜌0 𝑐02 exp −𝛼2 (𝑥2 + 𝑦2 ) (37)

Multiplying both sides of Eq. (49) by the smoothing function W (r 𝝃) yields ∫Ω

𝑊𝑖𝑗

⎧meshf ree method (SPH∕CSPM) ⎪ = ⎨meshf ree method (SPH∕CSPM) ⎪FDTD method ⎩

(29)

where 𝛼, 𝛽 = 1, 2, 3 are the dimension and 𝑓𝛼 (𝝃) =

𝑚𝑗 𝑗=1 𝜌𝑗

boundary t reat ment for particle 𝑖

The CSP method [36,37] improve the particle approximation in the standard SPH method using Taylor series. The Taylor series expansion of a function f (r) at a position 𝝃 can be represented as 1 𝛼 (𝒓 − 𝝃 𝛼 )(𝒓𝛽 − 𝝃 𝛽 )𝑓𝛼𝛽 (𝝃) + ... 2!

𝑓 (𝒓𝑗 )𝑊𝑖𝑗

∑𝑁

The finite-difference time-domain (FDTD) method is introduced to combine with the virtual particle technique, and thus a technique based on the meshfree-FDTD hybrid method for acoustic boundary treatment is accordingly constructed [48]. The feasibility and validity of the meshfree-FDTD hybrid method is verified by simulating sound propagation in one-dimensional pipes with boundaries. The present work use a two-dimensional version of the hybrid method. In the hybrid meshfree-FDTD boundary treatment, three types of particles need to be built before computation, namely the fluid particle, the boundary particle, and the virtual particle. Numerical method chosen for these three kinds of particles are shown as below

(28)

3.4. CSP method

𝑓 (𝒓) = 𝑓 (𝝃) + (𝒓𝛼 − 𝝃 𝛼 )𝑓𝛼 (𝝃) +

𝑚𝑗 𝑗=1 𝜌𝑗

4. Particle-based acoustic computation with different boundaries

𝑁 𝑁 ∑ ∑ 𝛿𝑝𝑖 + 𝛿𝑝𝑗 𝑚𝑗 𝐷0 𝛿𝒗𝑖 =− 𝑚𝑗 ∇𝑖 𝑊𝑖𝑗 + 𝛿𝒗𝑖 𝒗 ⋅ ∇𝑖 𝑊𝑖𝑗 𝐷𝑡 𝜌 𝜌 𝜌 𝑖𝑗 𝑖 𝑗 𝑗=1 𝑗=1 𝑗

+

(34)

𝑗

𝑁 𝑁 ∑ ∑ 𝑚𝑗 𝑚𝑗 𝐷0 𝛿𝜌𝑖 = 𝜌𝑖 𝛿𝒗𝑖𝑗 ⋅ ∇𝑖 𝑊𝑖𝑗 + 𝛿𝜌𝑖 𝒗 ⋅∇ 𝑊 𝐷𝑡 𝜌 𝜌 𝑖𝑗 𝑖 𝑖𝑗 𝑗=1 𝑗 𝑗=1 𝑗

− 𝛿𝒗𝑖 ⋅

.

In our Lagrangian approach, calculation of the particle motion and time integration is performed with a second-order leapfrog integration [47].

where 𝛿pi , 𝛿𝜌i , and 𝛿vi are the perturbations of pressure, density, and velocity of the fluid particle i due to acoustic waves. Considering the particle approximation of the gradient of the unity and the transformation method presented in [46], the particle approximation equations of the continuity and momentum equations of acoustic waves (Eqs. (26) and (27)) can be rewritten as

𝑁 ∑

∫Ω 𝑊 (𝒓−𝝃, ℎ)𝑑 𝒓

Following a similar derivation to that in the standard SPH method, the particle approximations for Eqs. (37) and (38) at particle i become

𝑁 ∑

𝑁 ∑ 𝑚𝑗

∫Ω 𝑓 (𝒓)𝑊 (𝒓−𝝃, ℎ)𝑑 𝒓

(𝒓𝛼 − 𝝃 𝛼 )𝑓𝛼 (𝝃)𝑊 (𝒓−𝝃, ℎ)𝑑 𝒓 462

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

Fig. 2. The sound pressure contour of sound propagation with rigid boundary at different time.

Fig. 3. Comparison of CSP results and exact solutions in solving sound propagation with the rigid boundary (a: sound pressure at x-axis; b: velocity perturbations at x-axis).

From the figure, it can be seen that, at y = 0 m, the CSP results possess similar smoothing trends compared to the theoretical solution. With the time increasing, the amplitude of the peaks and valleys decreases. Generally, the CSP algorithms and codes can simulate the Guassian pulse of sound propagation well. Sound propagation with perfectly matched layers (PML) is tested and results are shown in Fig. 4. Fig. 4(a) and (b) are sound pressure at t = 0.06 s and 0.09 s, respectively. Unlike rigid boundary, under the effect of the PML, the incident wave is absorbed with no reflected waves on the boundary. As shown in Fig. 5, CSP results are in good agreement with theoretical solutions, namely, the hybrid meshfree-FDTD method can accurately simulate the sound propagation with a PML boundary.

5.1. 1D sound propagation in steady uniform flow Sound propagation in steady uniform background flow is examined in this section. The test example is originally provided in the Fourth Computational Aeroacoustics Workshop on Benchmark Problems [49]. Consider sound propagation in a 1D homogeneous, steady uniform moving medium. The initial conditions are [ ] 𝛿𝑢(𝑥, 𝑡 = 0) = 𝛼1 𝑐0 [2 + cos(𝑘𝑥)] exp −𝛼2 (𝑘𝑥)2 , [ ] 𝛿𝑝(𝑥, 𝑡 = 0) = 𝛼1 𝜌0 𝑐02 [2 + cos(𝑘𝑥)] exp −𝛼2 (𝑘𝑥)2 , [ ] 𝛿𝜌(𝑥, 𝑡 = 0) = 𝛼1 𝜌0 [2 + cos(𝑘𝑥)] exp −𝛼2 (𝑘𝑥)2 ,

(38) (39) (40)

where 𝛿u is the velocity perturbation along the x-axis, t represents time, x is the geometric position, k is the wave number, and 𝛼 1 = 4/3402 , 𝛼 2 = ln 2/120. The fluid is an ideal gas with sound speed c0 340 m/s and density 1.0 kg/m3 . The exact solution for the linearized perturbation equations

5. Numerical tests and discussion In this section, different test cases with applications to modeling sound propagation in moving fluids are outlined. 463

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

Fig. 4. The sound pressure contour of sound propagation with the PML boundary at different time.

Fig. 5. Comparison of CSP results and exact solutions in solving sound propagation with PML boundary (a: sound pressure at x-axis; b: velocity perturbations at x-axis).

at time t is [ ] 𝛿𝑝(𝑥, 𝑡) = 𝛼1 𝜌0 𝑐02 [2 + cos(𝑘𝑥 − 𝑤𝑡)] exp −𝛼2 (𝑘𝑥 − 𝑤𝑡)2 ,

lines. The Mach number (Ma) of the background flow changes from −0.2 to 0.2. Results from the SPH method is not given here because they are so similar to the CSP method that the difference cannot be distinguished in the figure. The figure shows that sound pressure amplitude and tendency obtained by the CSP method agree well with exact solutions, and the CSP method is validated for this problem. To evaluate the convergence of the particle-based methods for computational acoustics, a non-dimensional sound pressure error (𝜀pre ) is used to discuss the effects of ∆x. 𝜀pre is defined as:

(41)

In the model, sound waves propagate along the x-axis over a time interval of 0.5 s, w is the angular frequency of 340 rad/s, and k is set as 1.0 rad/m. The computational domain spans −250 m to 250 m. Since the computational region is far longer than the sound propagation distance, Dirichlet boundary conditions are applied to both ends of the computational region. The sound pressure and perturbations of density and velocity for the last three particles at each end are set as 0, which are their initial value. In the computation, the particle spacing (∆x) is set to be 0.05 m, the Courant–Friedrichs–Lewy (CFL) number is 0.15, and the smoothing length (h) for the simulation is 1.45∆x. The CFL number is defined as Δt = CCFL Δx/c0 , where ∆x is the time step and CCFL is the CFL number. CSP results at time t = 0.5 s are compared with exact solutions in Fig. 6, and three test cases with the moving medium at different Mach numbers are shown. Computational results are plotted with points at intervals of 5 grid points, while exact solutions are plotted with solid

𝜀pre

√ √ 𝑁 √∑ ( )]2 √ [ ∗ 1 √ 𝑝𝑗 − 𝑝 𝑥𝑗 , = 2 3𝛼1 𝜌0 𝑐0 𝑁 𝑗=1

(42)

where N is the number of particles in the computational domain, 𝑝∗𝑗 is the numerical solutions of sound pressure for particle j, p(xj ) is the theoretical sound pressure at the position xj , and xj is the position of particle j. 464

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

Table 1 Comparison of SPH and CSP methods in solving PE and LPE at different Mach numbers. Mach number

Ma = −0.2 Ma = 0.0 Ma = +0.2

Equations

PE LPE PE LPE PE LPE

SPH

CSP

𝜀pre (× 10−5 )

tCPU (s)

𝜀pre (× 10−5 )

tCPU (s)

7.2831 7.6969 7.2828 7.6965 7.2824 7.6962

128 91 126 90 125 90

6.4546 6.8671 6.4543 6.8668 6.4539 6.8665

142 101 142 102 143 102

and the results are compared with theoretical solutions. The data show that 𝜀pre of the DNS method is lower than 1.0 × 10−8 , which is close enough to be considered exact. The bilinear interpolation is used to interpolate the DNS results to a value corresponding to the position of a fluid particle. As is shown in the table, both PE and LPE are validated to be useful in solving problems with acoustic propagation in steady uniform flow. Compared with PE, LPE can save some computational time with acceptable loss in the accuracy. Both two methods are successful in numerical implementation of the Lagrangian approach. Different Mach numbers of the uniform flow do little effects on 𝜀pre and tCPU .

Fig. 6. Sound pressure comparison between CSP and exact solutions for the LPE with the moving medium at different Mach numbers.

Convergence and CPU time (tCPU ) curves for SPH and CSP solutions of LPE are shown in Fig. 7. All computational parameters stay the same as in the model for Fig. 6 except ∆x. All computation in the present paper were performed with processor Intel(R) Core(TM) i7-4710MQ CPU 2.50 GHz and RAM 8.00 GB. According to our computations, the Mach number has little effects on the error, so only one Mach number is used in the following comparison. As shown in Fig. 7, two meshfree methods are convergent for this problem. The CSP method has better accuracy for different ∆x compared with the SPH method. As ∆x increases from 0.03 to 0.09 m, 𝜀pre for the SPH and CSP methods grow from about 2.0 × 10−5 to over 2.5 × 10−4 . The graph shows a fivefold increases for two methods in tCPU when ∆x decreases from 0.09 to 0.03 m. For further evaluation, Table 1 compares the accuracy and efficiency of the methods in solving perturbation equations and linearized perturbation equations. Since there is no theoretical solution for perturbation equations, the exact solutions come from DNS. In the DNS computation, we use a 4th order finite difference scheme for spatial derivative and a 4th order Runge–Kutta time integration scheme. In order to validate the DNS computation, it is used to solve linearized perturbation equations,

5.2. 1D sound propagation in unsteady flow After analyzing sound propagation in steady uniform flow, the same analysis is used to evaluate the Lagrangian approach in modeling sound propagation in unsteady flow. Theoretical solutions for a 1D unsteady background flow is given in Chapter 3 in [50]. Consider sound propagation in a 1D unsteady flow with 𝑢0 =

2 𝑥 + 𝐶1 , 𝛾 +1 𝑡

( )2 𝐶 2 (𝛾+1) 𝛾∕(𝛾−1) 2𝑥 1 ⎛ ⎞⎞ ⎛ 1 𝛾 − 1 ⎜− 2 𝐶1 + (𝛾+1)𝑡 + 2(𝛾−1) ⎟⎟ ⎜ 𝑝0 = 𝐶3 , 2−2𝛾 ⎜ 𝐶3 𝛾 ⎜ 𝐶2 (𝛾−3) 𝛾+1 ⎟⎟ 𝑥2 + (𝛾+1 ⎝ ⎝+ 𝛾+1 𝑡 ⎠⎠ )𝑡2

(44)

( )2 𝐶 2 (𝛾+1) 1∕(𝛾−1) 2𝑥 ⎛ 1 ⎛ ⎞⎞ (𝛾 − 1) ⎜− 2 𝐶1 + (𝛾+1)𝑡 + 21(𝛾−1) ⎟⎟ ⎜ 𝜌0 = , ⎜ 𝐶3 𝛾 ⎜+ 𝐶2 (𝛾−3) 𝑡(2−2𝛾)∕(𝛾+1) + 𝑥2 ⎟⎟ ⎝ ⎝ 𝛾+1 (𝛾+1)𝑡2 ⎠⎠

(45)

Fig. 7. Convergence and CPU time curves for SPH and CSP solutions of LPE at Ma = 0.2. 465

(43)

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

Fig. 8. Pressure of CSP solutions of LPE versus exact solutions. (a: pressure at different time (perturbations are sound pressure); b: detailed view of dashed box in (a)).

Fig. 9. Convergence and CPU time curves for SPH and CSP simulation of sound propagation in unsteady flow.

( 𝑐02 =

𝛾 −1 𝑥 − 𝐶1 𝛾 +1 𝑡

)2 −

3−𝛾 (𝛾 − 1)𝐶2 𝑡−2(𝛾−1)∕(𝛾+1) , 𝛾 +1

For further evaluation, convergence and CPU time curves for SPH and CSP solutions of LPE are shown in Fig. 9. All computational parameters stay the same as in the model for Fig. 8 except ∆x. Fig. 9 shows that the SPH method does not converge, while the CSP method does. This results consistent with the analysis in [51], which gives an error estimation for four schemes including both the standard SPH method and the CSP method. For irregular particle arrangement, [51] shows that the CSP method converges as ∆x decreases, while convergence for the standard SPH method requires both ∆x → 0 and h / ∆x → ∞. This conclusion can also be found in [52]. On the other hand, for the CSP method, 𝜀pre decreases along with ∆x. A decrease of ∆x increases the time cost as can be seen from Fig. 9. tCPU increases from less than 100 s to about 600 s when ∆x decreases from 0.09 m to 0.03 m. The SPH method has the lowest tCPU cost. Nevertheless, the difference is small between the two methods. For further evaluation, Table 2 compares 𝜀pre and tCPU of the SPH and CSP methods in solving PE and LPE with unsteady background flow. The table indicates that both PE and LPE are able to solve acoustic problems with unsteady background flow. In addition, LPE save much computational time with acceptable loss in the accuracy comparing with PE. In solving PE, similarly to the previous results of solving LPE in Fig. 7, the SPH method cost the least computational time.

(46)

where u0 is the initial velocity of background flow, C1 = 30.0, C2 = −1.0 × 106 , C3 = 82,571.0, 𝛾 =1.4, and the initial time (t0 ) is at 12.5 s. By selecting these values, the initial velocity of background flow is about 40 m/s, and the sound speed is around 340 m/s. Both depend on position. Sound waves propagate along the x-axis over a time interval of 0.5 s, and the initial sound waves are set as Eqs. (38–40). The computational domain spans −250 m to 250 m. The same Dirichlet boundary conditions to Section 5.1 are applied to both ends of the computational region. In the computation, Δx is set to be 0.05 m, CCFL is 0.15, and h is selected as 1.45Δx. A comparison of pressure obtained from CSP results in solving LPE and exact solutions is given in Fig. 8, and Fig. 8(b) gives a detailed view of the dashed box in Fig. 8(a). The CSP results are plotted with points at intervals of 30 grid points in figure (a) and 15 grid points in figure (b), and the exact solutions are plotted with solid lines. The exact solutions come from DNS results (details see Section 5.1). The graph illustrates that CSP results get good agreement with exact solutions in solving LPE. 466

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

Fig. 10. Contour of density change obtained by the CSP simulation of sound propagation in moving fluid at different time.

Fig. 11. Comparison between CSP results and exact solutions along the x-axis at 0.05 s (a: sound pressure; b: velocity perturbations). Table 2 Comparison of SPH and CSP methods in solving PE and LPE with unsteady flow. Equations

SPH

CSP

𝜀pre (× 10 )

tCPU (s)

𝜀pre (× 10−5 )

tCPU (s)

27.1432 27.5736

267 195

5.8515 6.3053

284 207

−5

PE LPE

computational aeroacoustics (CAA) [53]. In addition, the standard SPH method is not discussed in 2D problems, because it has more limitation to maintain convergence and the selection of h is also more difficult than the CSP method in the application as shown in Sections 5.1 and 5.2. Consider sound propagation in steady uniform √ flow with Ma = 0.2. Let 𝛼 1 = 10/3402 , 𝛼 2 = ln 2/9, 𝛼 3 = ln 2/7, 𝜂 = (𝑥2 + 𝑦2 ), where x and y are the position coordinates. The initial conditions are [ ( )] 𝛿𝑝 = 𝛼1 𝜌0 𝑐02 exp −𝛼2 𝑥2 + 𝑦2 , [ ( 2 )] [ ( )] 𝛿𝜌 = 𝛼1 exp −𝛼2 𝑥 + 𝑦2 + 0.1𝛼1 exp −𝛼3 (𝑥 − 10)2 + 𝑦2 , [ ( )] 𝛿𝑢 = 0.04𝛼1 𝑐0 𝑦 exp −𝛼3 (𝑥 − 10)2 + 𝑦2 , [ ( )] 𝛿𝑣 = −0.04𝛼1 𝑐0 (𝑥 − 10) exp −𝛼3 (𝑥 − 10)2 + 𝑦2 ,

5.3. 2D sound propagation in steady uniform flow In this case, a 2D acoustic model is selected to validate the Lagrangian approach. This model is similar to the problem given in the ICASE/LaRC workshop concerning benchmark problems in 467

(47) (48) (49) (50)

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

Fig. 12. Convergence curves for CSP solutions of LPE at Ma = 0.2.

Fig. 13. Contour of CSP solutions of density change at different time.

where 𝛿 u and 𝛿 v are the velocity perturbation of particles along the xaxis and y-axis, respectively. The fluid is an ideal gas with sound speed c0 340 m/s and 𝜌0 1.0 kg/m3 . The exact solutions for the sound waves when time equals t are ( 2) ∞ 𝛼1 𝜌0 𝑐02 −𝜉 𝛿𝑝 = exp cos (𝜉𝑡)𝐽0 (𝜉𝜂)𝜉𝑑𝜉, (51) 2𝛼2 ∫0 4 𝛼2

taken to be Dirichlet boundary with perturbation variables set to 0. A model with larger computational domain is implemented to show the process of sound propagation in the uniform flow at first, then a smaller model is used for the evaluation in order to save the time. The CSP results for the change in density at different time steps are shown in Fig. 10. In this simulation, the computational domain is −80 m < x < 80 m, −80 m < y < 80 m. During the computation, all computational particles are moving along with the background flow, and sound waves propagate in the media at the same time. The direction of the background flow is represented by an arrow. Computational

where J0 () and J1 () are Bessel functions of order 0 and 1, respectively. In the simulation, ∆x is set as 0.2 m, CCFL is selected as 0.2, and h is 1.1∆x. The boundary particles around the computational domain are

468

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

Fig. 14. Comparison of density change between CSP results and exact solutions along y = 10 m.

Fig. 15. Convergence and CPU time curves for CSP solutions of sound propagation in vortex flow.

results show that the CSP method can simulate the propagation process of Gaussian impulses. To evaluate the computational results, Fig. 11 gives a comparison between CSP results and DNS solutions along the x-axis. In this simulation, the computational domain is −40 m < x < 40 m, −40 m < y < 40 m, and the time for sound propagation is 50 ms. The points in the figure represent the numerical values at different positions, and the lines denote the exact solutions. In order to clearly identify numerical results, the points are plotted at intervals of 5 grid points. As can be seen from the line graph, numerical results of all perturbation variables agrees well with exact solutions. The comparison suggests that the CSP results are accurate. Convergence and CPU time curves for CSP solutions of LPE are shown in Fig. 12. As can be seen from the figure, the CSP method converge when solving LPE. 𝜀pre increases from about 2.0 × 10−6 to more than 1.0 × 10−4 along with ∆x increasing from 0.1 to 0.4 m. At the same time, tCPU decreases significantly along as ∆x increase. It changes from more than 1000 s to around 20 s when ∆x changes from 0.10 m to 0.40 m. Numerical results show that both PE and LPE are able to solve 2D acoustic problems in steady uniform flow. 𝜀pre and CPU time for PE and LPE solving with CSP are 1.6742 and 1.6805, 221 s and 130 s

respectively. LPE can save some computational time with acceptable loss in the accuracy comparing with PE. 5.4. 2D sound propagation in rotational vortex 2D sound propagation in a rotational vortex is considered in this case. The vortex is described by u0 = −𝜔y, v0 = 𝜔x, 𝑝0 = 12 𝜔2 𝑥2 + 1 2 2 𝜔 𝑦 2

+ 𝐶pre , where 𝜔 is the angular velocity of the rotational vortex, x and y are the position coordinates, and u0 and v0 are the velocity of particles along the x-axis and y-axis, respectively. A 2D acoustic model similar to Section 5.3 is selected to validate the Lagrangian approach. Let 𝛼1 = 102 , 𝛼2 = ln92 , 𝛼3 = ln72 . 340 The initial conditions are { [ ]} 𝛿𝑝 = 𝛼1 𝜌0 𝑐02 exp −𝛼2 (𝑥 − 10)2 + (𝑦 − 10)2 , (52) { [ ]} 𝛿𝜌 = 𝛼1 exp −𝛼2 (𝑥 − 10)2 + (𝑦 − 10)2 { [ ]} + 0.1𝛼1 exp −𝛼3 (𝑥 − 20)2 + (𝑦 − 10)2 ,

(53)

A model with larger computational domain is implemented to show the process of sound propagation in the rotational flow at first, then a smaller model is used for the evaluation in order to save the time. 469

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

The contour of density change at different time is shown in Fig. 13. In this simulation, the computational domain is −30 m < x < 130 m, −30 m < y < 130 m, and 𝜔 is set as 8 rad/s. As shown in the contour, the red line at position (20, 10) is the center of the rotational vortex. All particles are moving during the computation, and sound waves propagate in the media at the same time. For further evaluation, a model close to the actual situation is considered. The computational domain is −30 m < x < 50 m, −30 m < y < 50 m, and 𝜔 is set as 2 rad/s. Fig. 14 gives a comparison of density change between the CSP method and exact solutions at y = 10 m. The line represents exact solutions, and points in the figure represent the CSP results. As can be seen from the line graph, the CSP results show good agreement with the exact solutions. Convergence and CPU time curves for CSP solutions of LPE are shown in Fig. 15. All computational parameters stay the same as in Section 5.3 except ∆x. Fig. 15 shows that the CSP method are convergent in this case. 𝜀pre for two meshfree methods increases from less than 1.0 × 10−5 to over 1.0 × 10−4 along with ∆x changing from 0.1 to 0.4 m. For the CSP method, t CPU increases over 10 times when ∆x decreases from 0.4 to 0.1 m. For further evaluation, 𝜀pre and tCPU of CSP method in solving PE and LPE are 1.8549 and 1.8612, 223 and 136 respectively when ∆x = 0.2 m.

[7] Li E, He ZC. Development of a perfect match system in the improvement of eigenfrequencies of free vibration. Appl Math Model 2017;44:614–39. [8] Chai YB, Gong ZX, Li W, Li TY, Zhang QF, Zou ZH, Sun YB. Application of smoothed finite element method to two-dimensional exterior problems of acoustic radiation. Int J Comput Methods 2018;15:1850029. [9] Chai YB, You XY, Li W, Huang Y, Yue ZJ, Wang MS. Application of the edge-based gradient smoothing technique to acoustic radiation and acoustic scattering from rigid and elastic structures in two dimensions. Comput Struct 2018;203: 43–58. [10] Uras RA, Chang CT, Chen Y, Liu WK. Multiresolution reproducing kernel particle methods in acoustics. J Comput Acoust 1997;5:71–94. [11] Bouillard P, Suleaub S. Element-free Galerkin solutions for Helmholtz problems: fomulation and numerical assessment of the pollution effect. Comput Method Appl Mech Eng 1998;162:317–35. [12] Karageorghis A, Lesnic D, Marin L. The MFS for the identification of a sound-soft interior acoustic scatterer. Eng Anal Bound Elem 2017;83:107–12. [13] He Z, Li P, Zhao G, Chen H. A meshless Galerkin least-square method for the Helmholtz equation. Eng Anal Bound Elem 2011;35:868–78. [14] Fu Z, Chen W, Gu Y. Burton–Miller-type singular boundary method for acoustic radiation and scattering. J Sound Vib 2014;333:3776–93. [15] Li J, Fu Z, Chen W. Numerical investigation on the obliquely incident water wave passing through the submerged breakwater by singular boundary method. Comput Math Appl 2016;71:381–90. [16] Fu Z, Chen W, Wen P, Zhang C. Singular boundary method for wave propagation analysis in periodic structures. J Sound Vib 2018;425:170–88. [17] Tang Z, Fu Z, Zheng D, Huang J. Singular boundary method to simulate scattering of SH wave by the canyon topography. Adv Appl Math Mech 2018;10:912–24. [18] Leighton T. The acoustic bubble. London: Academic Press; 2012. [19] Schwarz A, Janicka J. Combustion noise, 102. Berlin: Springer Science & Business Media; 2009. [20] Kytömaa HK. Theory of sound propagation in suspensions: a guide to particle size and concentration characterization. Powder Technol 1995;82:115–21. [21] Liu MB, Liu GR. Smoothed particle hydrodynamics (SPH): an overview and recent developments. Arch Comput Methods Eng 2010;17:25–76. [22] Liu MB, Liu GR, Zong Z. An overview on smoothed particle hydrodynamics. Int J Comput Methods 2008;5:135–88. [23] Monaghan JJ. Smoothed particle hydrodynamics and its diverse applications. Annu Rev Fluid Mech 2012;44:323–46. [24] Wolfe CT, Semwal SK. Acoustic modeling of reverberation using smoothed particle hydrodynamics Master Thesis. University of Colorado; 2008. [25] Hahn P, Negrut D. On the use of meshless methods in acoustic simulations. In: Proceedings of the ASME 2009 international mechanical engineering congress and exposition, IMECE2009; 2009. p. 185–99. [26] Treyssede F, Gabard G, Tahar MB. A mixed finite element method for acoustic wave propagation in moving fluids based on an Eulerian–Lagrangian description. J Acoust Soc Am 2003;113:705–16. [27] Gabard G, Lefrancois E, Tahar MB. Aeroacoustic noise source simulations based on Galbrun’s equation. In: Proceedings of the 10th AIAA/CEAS aeroacoustics conference, Manchester, Great Britain. AIAA; 2004. 2004-2892. [28] Brekhovskikh LM, Godin OA. Acoustics of layered media II: point sources and bounded beams, 2. Berlin, Heidelberg: Springer; 1999. [29] Seo JH, Moon YJ. Perturbed compressible equations for aeroacoustic noise prediction at low mach numbers. AIAA J 2005;43:1716–24. [30] Seo JH, Moon YJ. Linearized perturbed compressible equations for low Mach number aeroacoustics. J Comput Phys 2006;218:702–19. [31] Ewert R, Schröder W. Acoustic perturbation equations based on flow decomposition via source filtering. J Comput Phys 2003;188:365–98. [32] Lucy LB. A numerical approach to the testing of the fission hypothesis. Astron J 1977;82:1013–24. [33] Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon Not R Astron Soc 1977;181:375–89. [34] Springel V. Smoothed particle hydrodynamics in astrophysics. Annu Rev Astron Astrophys 2010;48:391–430. [35] Liu GR, Liu MB. Smoothed particle hydrodynamics: a meshfree particle method. Singapore: World Scientific; 2003. [36] Chen JK, Beraun JE, Carney TC. A corrective smoothed particle method for boundary value problems in heat conduction. Comput Methods Appl Mech Eng 1999;46:231–52. [37] Chen JK, Beraun JE, Jih CJ. Completeness of corrective smoothed particle method for linear elastodynamics. Comput Mech 1999;24:273–85. [38] Zhang ZL, Liu MB. A decoupled finite particle method for modeling incompressible flows with free surfaces. Appl Math Model 2018;60:606–33. [39] Liu MB, Xie WP, Liu GR. Modeling incompressible flows using a finite particle method. Appl Math Model 2005;29:1252–70. [40] Koshizuka S, Shibata K, Kondo M, Matsunaga T. Moving particle semi-implicit method: a meshfree particle method for fluid dynamics. Cambridge: Academic Press; 2018. [41] Colagrossi A, Landrini M. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. J Comput Phys 2003;191:448–75. [42] Monaghan JJ. Simulating free surface flows with SPH. J Comput Phys 1994;110:399–406. [43] Monaghan JJ. Implicit SPH drag and dusty gas dynamics. J Comput Phys 1997;138:801–20. [44] Zhang YO, Zhang T, Ouyang H, Li TY. SPH simulation of sound propagation and interference. In: Proceedings of the 5th international conference on computational method, ICCM’14, Cambridge; 2014.

6. Conclusions A Lagrangian approach for computational acoustics is presented and implemented numerically, using the SPH and CSP methods to solve perturbation equations. In this Lagrangian approach, both acoustic and fluid variables are supposed to be the physical properties of a fluid particle which can move through both space and time, and the perturbation equations are derived based on the hydrodynamic/acoustic splitting method. Applications to modeling sound propagation in moving fluids are outlined with several examples, and a comparison of the PE and LPE using the particle-based methods is presented. Effects of particle spacing on SPH and CSP methods are evaluated. The error analysis show that the present Lagrangian approach gives good agreement with exact solutions, and LPE can save some computational time with acceptable accuracy loss compared with PE. For different test cases, the particle-based methods are accurate for the time domain simulation of acoustic waves. The CSP method is convergent for all cases considered, while the standard SPH method is not. As illustrated in this paper, the Lagrangian approach is supposed to have potential for acoustic problems with moving boundaries and multiphase systems like bubble acoustics. Nevertheless, before engineering applications, variety of tests for the present Lagrangian approach still need to be analyzed. This issue is still under study and will be reported in a subsequent paper. Acknowledgments This work was supported by the National Natural Science Foundation of China (Nos. 51809208 and 51741910). References [1] Wang SZ. Finite-difference time-domain approach to underwater acoustic scattering problems. J Acoust Soc Am 1996;99:1924–31. [2] Harari I. A survey of finite element methods for time-harmonic acoustics. Comput Methods Appl Mech Eng 2006;195:1594–607. [3] Kythe PK. An introduction to boundary element methods, 4. Boca Raton: CRC Press; 1995. [4] He ZC, Li GY, Zhong ZH, Cheng AG, Zhang GY, Li E, Liu GR. An ES-FEM for accurate analysis of 3D mid-frequency acoustics using tetrahedron mesh. Comput Struct 2012;106:125–34. [5] He ZC, Li E, Li GY, Wu F, Liu GR, Nie X. Acoustic simulation using 𝛼-FEM with a general approach for reducing dispersion error. Eng Anal Bound Elem 2015;61:241–53. [6] Li E, He ZC, Xu X, Liu GR. Hybrid smoothed finite element method for acoustic problems. Comput Method Appl Mech Eng 2015;283:664–88. 470

Y. Zhang, S.G.L. Smith and T. Zhang et al.

Engineering Analysis with Boundary Elements 108 (2019) 459–471

[45] Zhang YO, Zhang T, Ouyang H, Li TY. SPH simulation of acoustic waves: effects of frequency, sound pressure, and particle spacing. Math Probl Eng 2015;7:348314. [46] Zhang YO, Zhang T, Ouyang H, Li TY. Efficient SPH simulation of time-domain acoustic wave propagation. Eng Anal Bound Elem 2016;62:112–22. [47] Kelager M. Lagrangian fluid dynamics using smoothed particle hydrodynamics Ph.D. Thesis. University of Copenhagen; 2006. [48] Zhang YO, Li X, Zhang T. Modeling sound propagation using the corrective smoothed particle method with an acoustic boundary treatment technique. Math Comput Appl 2017;22:26.

[49] Dahl M.D. Fourth computational aeroacoustics (CAA) workshop on benchmark problems. NASA/CP-2004-212954, E-14393; 2004. [50] Von Mises R. Mathematical theory of compressible fluid flow. New York: Dover Publications; 2004. [51] Fatehi R, Manzari MT. Error estimation in smoothed particle hydrodynamics and a new scheme for second derivatives. Comput Math Appl 2011;61:482–98. [52] Zhu Q, Hernquist L, Li Y. Numerical convergence in smoothed particle hydrodynamics. Astrophys J 2015;800:6. [53] Hardin JC, Ristorcelli JR, Tam CKW. ICASE/LaRC workshop on benchmark problems in computational aeroacoustics. NASA CP-3300, 1995.

471