Simulation of the acoustic wave propagation using a meshless method

This paper presents numerical simulations of the acoustic wave propagation phenomenon modelled via Linearized Euler equations. A meshless method based on collocation of the strong form of the equation system is adopted. Moreover, the Weighted least squares method is used for local approximation of derivatives as well as stabilization technique in a form of spatial filtering. The accuracy and robustness of the method is examined on several benchmark problems.


Introduction
Numerical study of sound generation and propagation phenomenon plays a major role in Computational Aeroacoustics (CAA), especially in predictions of sound generated by turbulent flows and its environmental consequences.This paper is aiming to contribute to the development of suitable methods for sound propagation simulations governed by the Linearized Euler equations (LEE), cf.[1,2].
Various numerical mesh-based methods, cf.[3,4], have been investigated in this context, but there is also growing interest in meshless methods, cf.[5][6][7][8][9][10][11], due to their potential advantages.Depending on applications, these methods can be superior in accuracy, stability, robustness and can deal with complex geometries, multiphase flow, problems with free boundaries and moving objects.In this paper, the development of a meshless method adjusted for linear systems is outlined.
Firstly, LEE as a linear hyperbolic system is introduced, followed by the Weighted Least SQuares (WLSQ) method which is employed in two ways -for local interpolation of acoustic variables and for spatial filtering, cf.[12,13].The WLSQ filtration technique serves primarily for the stabilization of the numerical scheme and can be utilized under the assumption that the governing equations are linear and the solution is smooth enough.
Finally, abilities of the method are demonstrated on various benchmark problems for the acoustic wave propagation, cf.[14][15][16].The method is tested on unstructured point distributions and several types of boundary conditions (BC) are implemented, namely the nonreflecting BC realized by the Perfectly Matched Layer (PML), cf.[17][18][19].

Linearized Euler equations
The hyperbolic non-homogeneous LEE represent one of the suitable models for the acoustic wave propagation phenomenon.We consider the 2D matrix form (x = (x, y) ∈ Ω, t > 0) written as where Vector function w (x, t) denotes the time dependent acoustic variables (density, velocity components and pressure) and w 0 (x) steady mean flow variables corresponding to the underlying flow field.Jacobian matrices of the system (1) are given as where γ = 1.4 is the ratio of specific heats and S = S(x, t) represents the acoustic source term, cf.[1,2].The initial-value boundary problem for LEE then consists of the equation system (1) with the initial and boundary condition

Meshless method with WLSQ filtering
In the context of meshless methods we distiguish between the global cloud Ω defined as a set of n points discretizing a domain of interest Ω and local clouds Ωi , i = 1, . . ., n that are calculated in the pre-processing step of the simulation.An i-th local cloud usually consists of a set of n i neighbours chosen by a certain rule.Practically, efficient algorithms for nearest neighbor search such as the k-d tree are used.Therefore, the number of points n i is controlled by the radius of an open ball or it can be prescribed manually.

Weighted Least Squares Method
For every local cloud Ωi we wish to find a local approximation ŵ : where and Functions p l : R d → R, l = 1, . . ., m, form a basis B of an approximation space F = span(B).The complete (multivariate) polynomial basis of degree ν, cf.[5,10] is adopted in this work.For d = 2 and ν = 3, the basis

WLSQ approximation
It is well known that WLSQ method finds the coefficients α in ( 6) by minimizing the objective function J(α) with weights where P denotes the moment matrix of type (n i × m), Φ the diagonal weight matrix of order n i and w the vector of given function values.The weighting function φ : R d → R prescribes the contribution of every point to J(α) in the local cloud according the distance to the star point x i 1 , cf. [10].Solution to the system of normal equations P T ΦPα = P T Φw. (11) then yields α = P T ΦP −1 P T Φw (12) assumming the invertibility of the above mentioned matrix.We suppose that the number of points n i in Ωi is greater than the number of basis functions m, i.e. n i > m holds.

WLSQ interpolation
By adding the interpolation constraint to the Normal equations (11), the value w 1 is reproduced at the star point x i exactly.This modification can be gently implemented and allows to use the WLSQ method for multiple purposes, namely the spatial discretization of the governing equations (interpolation) as well as numerical filtering (approximation).

Spatial discretization of LEE
In order to obtain the semi-discrete form of governing equations (1), the collocation method is adopted, i.e. the expression has to be satisfied.Acoustic variables ŵi are locally approximated using WLSQ method with interpolation condition described in section (3.1).
The collocation then leads to a system of ordinary differential equations (ODE) in form where RHS i (t) contains contributions to the righthand side from all neighbours of the star point x i , cf. [7,10,11].The interpolation condition in the WLSQ approach allows to simplify the left side of the equation system (15) avoiding the necessity to solve the system of n linear equations.Therefore, using WLSQ with interpolation condition greatly reduces the computation costs and improves stability, cf.[10,11].The system of ODE ( 15) is then solved with high order lowdissipation and low-dispersion runge-kutta scheme optimized for wave propagation problems, cf.[20].

Numerical filtering
The numerical scheme (15) can be viewed as a generalized finite differences scheme based on central differences which is not stable in general, cf.[8,9].For linear problems, e.g. the acoustic wave propagation described by LEE, the stabilization can be gently achieved in the context of WLSQ method used as a filter.In general, filters based on least squares minimization are known as Savitzky-Golay filters, cf.[12,13].It is not only the stability, but also robustness and accuracy of the meshless method is then affected by the WLSQ filter.
The i-th local cloud Ωi is utilized for filtering after every time step in a way that w i (t) ≈ w (x i , t) is replaced by the filtered value w i (t) = p T (x i )α ( 16) at the star point x i .

Boundary conditions
Firstly, solid-wall slip BC is prescribed at points where the acoustic waves should be reflected.The acoustic velocity vector v = (u , v ) at the boundary point x ∈ Γ is forced to stay orthogonal to the normal vector n = (n x , n y ), i.e.
Secondly, non-reflecting BC is prescribed at domain boundaries, where waves should leave the computational domain without any reflection.It can be achieved artificially by designing a layer of points around the domain of interest that will affect the solution in such a way that the major spectrum of incoming waves will dissipate.The PML is implemented in our numerical experiments, cf.[17][18][19].

Numerical experiments 4.1 2D Wall-reflected acoustic pulse problem
This benchmark problem is designed to verify the implementation of the wall boundary condition and the accuracy and stability of a numerical method when a reflection of an acoustic wave from the solid wall occur.Let us consider a 2D domain Ω = (−100, 100) × 0, 150) that is depicted in Fig. 1 together with an initial pulse located at (x p , y p ) = (0, 25).The wall boundary condition is prescribed at line y = 0. We suppose the underlying medium at rest, which implies that the mean flow velocity equals zero and therefore the Mach number M = 0 and for w 0 it holds The Initial-value boundary problem consists of the equation system (1) with the source term S = 0, above mentioned boundary conditions and initialized by the initial condition w (x, y, 0) given as follows where the radius r p = (x − x p ) 2 + (y − y p ) 2 and κ = (ln 2)/b 2 .The amplitude and the half-width of the acoustic initial pulse are determined by ε = 1 and b = 5, respectively.The comparison of the acoustic pressure p (x, y, T ) along the line x = 0 and analytical solution at time T = 75 is shown in Fig. 6.

2D Convected acoustic monopole
The domain of interest for the following benchmark problem is simple square Ω = (−100, 100) 2 extended by the PML Ω e , cf.Fig. 7. Furthemore, stretching of points towards outlet boundaries and artificial dissipation is applied in order to suppress all incoming waves and eliminate spurious reflections.We will consider subsonic underlying mean flow in x-direction with Mach number M 0 = 0.5 as well as supersonic flow with M 0 = 1.5.In this example, the acoustic source term on the right-hand side of the system (1) is the periodic perturbation of acoustic density and pressure at the origin x s = y s = 0 which is given as S(x, t) = e −α((x−xs) 2 +(y−ys

2D Acoustic monopole in a sheared mean flow
Prescribing the velocity components in (18) in form we get the non-constant underlying flow known as sheared mean flow, cf.[1,15].In this example, the peak velocity is taken as u max = 0.5 and the shear layer thickness δ = 150.The convection velocity profile is depicted in Fig. 10.Instantaneous acoustic pressure contours at time T = 270 are plotted in Fig. 11.

Discussion
Linearity of governing equations (LEE) and therefore the absense of shock waves and strong disturbances in the solution allowed to develop a new stabilization technique in a meshless framework.Instead of the upwind approximation of derivatives widely used, e.g. in Finite volume methods, the stability of the meshless method was achieved by the WLSQ filtering of the numerical solution after each time step.
The method was tested on acoustic benchmark problems, such as the 2D wall bounded acoustic pulse problem and the acoustic monopole in a free stream and sheared mean flow.In order to suppress waves leaving the domain, the non-reflecting boundary condition in form of the PML was sucessfully utilized, cf.Fig. 8.To conclude, the method proved its ability to solve the wave propagation problems governed by LEE in various settings and is prepared for further applications.

4. 1 . 1
Numerical solution Results of the acoustic wave simulation in terms of the acoustic pressure contours are plotted ase a series of Fig. (2-5) at particular times T = 10, 25, 50, 75.The time integration was performed with time step dt = 0.025 that corresponds to the unstructured point distribution with n = 36006 number of points.
) 2 ) sin(ωt)(1, 0, 0, 1) T (20) where α = ln(2)/2, amplitude ε = 0.5 and angular frequency ω = 2π/30.The temporal computation was performed with time step dt = 0.1, number of points n = 50181 irregularly distributed and stretched in the PML.Contours of the acoustic pressure p (x, y, T ) at time T = 270 for the subsonic case are shown in Fig. 8 and for the supersonic case in Fig. 9. Waves moving in the PML layer clearly dissipate and there are no spurious waves returning to the domain of interest.