Analysis of Aeroacoustic Properties of the Local Radial Point Interpolation Cumulant Lattice Boltzmann Method

The lattice Boltzmann method (LBM) has recently been used to simulate wave propagation, one of the challenging aspects of wind turbine modeling and simulation. However, standard LB methods suffer from the instability that occurs at low viscosities and from its characteristic lattice uniformity, which results in issues of accuracy and computational efficiency following mesh refinement. The local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM) is proposed in this paper to overcome these shortcomings. The LB equation is divided into collision and streaming steps. The collision step is modeled by the cumulant method, one of the stable LB methods at low viscosities. In addition, the streaming step, which is naturally a pure advection equation, is discretized in time and space using the Lax–Wendroff scheme and the local radial point interpolation method (RPIM), a mesh free method. We describe the propagation of planar acoustic waves, including the temporal decay of a standing plane wave and the spatial decay of a planar acoustic pulse. The analysis of these specific benchmark problems has yielded qualitative and quantitative data on acoustic dispersion and dissipation, and their deviation from analytical results demonstrates the accuracy of the method. We found that the LRPIC-LBM replicates the analytical results for different viscosities, and the errors of the fundamental acoustic properties are negligible, even for quite low resolutions. Thus, this method may constitute a useful platform for effectively predicting complex engineering problems such as wind turbine simulations, without parameter dependencies such as the number of points per wavelength Nppw and resolution σ or the detrimental effect caused by the use of coarse grids found in other accurate and stable LB models.


Introduction
Evidence of early sailing boats on the Nile and of Persian pumps and mills from the first century B.C. shows humans have been interested in Wind Energy since ancient times [1]. In general, a wind turbine is defined as a device which converts the wind's kinetic energy into electrical energy [2,3]. It plays a key role on producing intermittent renewable energy and implementing a strategy to lower costs and reducing the reliance on fossil fuels [4,5]. Wind turbines have unique aerodynamic and aeroacoustic behavior that makes their prediction most challenging [6,7], particularly their simulation needs an enormous number of grid points or cells, and long enough time samples [8]. Researchers and centers such as the National Renewable Energy Laboratory (NREL) and the National Wind Technology Center (NWTC) have initiated multi-year programs on aeroacoustic wind turbine modeling [9] to develop efficient and appropriate computational aeroacoustic (CAA) implementations. Among particular issues specified to wind turbine problems, the propagation of sound is always a significant computational challenge [10,11]. With this aim, different numerical approaches were developed in the field of computational aeroacoustics. Tam [12] and Wells et al. [13] proposed popular numerical schemes such as compact and non-compact optimized schemes like the high-order compact difference scheme [14] and the dispersion-relation-preserving (DRP) scheme [15]. Cheong et al. proposed grid-optimized dispersion-relation-preserving (GODRP) schemes for aeroacoustic simulations [16]. Due to the huge cost of CAA simulations, hybrid methods using two sets of equations, one for the flow and another one for the acoustic disturbance field were developed [17].
Direct aerostatic simulations are computer-intensive due to the small ratio between sound pressure and pressure variation as a whole, the spreading of acoustic fields over a large area, and the time-consuming nature of traditional methods [18]. As an illustration, the direct numerical simulation of waves using Navier-Stokes equations requires schemes of fifth-order accuracy in space and fourth order accuracy in time [19,20]. Therefore, the lattice Boltzmann method (LBM), an explicit time marching scheme [21], has been widely used as an alternative to simulate sound wave propagation due to its kinetic nature, relative simplicity of implementation and parallelization. The LBM as a mesoscopic method uses probability density functions (probability of finding particles within a certain range of velocities at a certain range of locations at a given time) to model the momentum distribution in discrete space, thereby economizing computer resources [22]. Buick et al. [23] and Dellar et al. [24] studied sound wave propagation using LBM and achieved acceptable results. Bres et al. [25] and Gorakifard et al. [26] presented the dissipation and dispersion of acoustic waves using the BGK-LBM and the cumulant LBM, respectively. Furthermore, a regularized method for the BGK-LBM [27] and the recursive and regularized LBM (LBM-rrBGK) [28,29] have been developed to model wave propagation.
One of the drawbacks of the LBM is lattice uniformity, originated from symmetric lattice velocity models such as the square and cubic lattice meshes in 2D/3D simulations [30]. Lattice uniformity causes the streaming step to occur at uniform neighboring grid points. Thus, the LBM is only applied to uniform meshes, whereas issues of accuracy and computational efficiency mainly affect simulations of problems that require non-uniform meshes. For example, numerical simulations with curved and irregular boundaries, common in wind turbine simulations, encounter difficulties when fitting the grids to the boundaries or adapting complex computational domains. Grid refinement schemes or adaptive LBM can help to simulate curved and moving boundaries more accurately [31]. Wood [32] used refinement in LBM simulations to analyze wind energy and utilized adaptive LBM for moving boundaries [33]. However, these schemes are associated with higher computational costs and even additional perturbations in acoustical problems [34].
Generally, the methods commonly used along with the LBM on non-uniform meshes include three distinct categories [35]. The first is the interpolation-supplemented LBM (IS-LBM) [36,37]. This method adds an interpolation step to the collision and streaming steps of the LBM. Two major drawbacks of the IS-LBM are its inefficiency due to the need to interpolate at each time-step, and the appearance of negative particle distributions [38]. The second method is the combination of the LBM with the finite difference method (FDLBE) [39], finite volume methods (FVLBE) [40][41][42], or finite element methods (FELBE) [43][44][45] used to stabilize the computation. The third method is the Taylor-series expansion and least-squares-based lattice Boltzmann method (TLLBM) [46,47] instead of direct interpolation. These methods, the implementation of which is somewhat simple, use continuous distribution functions in physical space.
Although the above numerical schemes show robustness for complex problems, they are affected by the inherent shortcomings of using meshes in numerical methods, such as the enormous cost of generating meshes, the low level of stress accuracy in fluid-structure interaction simulations (FSI) [48], obstacles in the adaptive analysis, and limitations in simulations of physical phenomena with singular, or nearly singular behavior. Mesh-free or meshless methods have been devised to eliminate mesh-related problems [49]. The MFree method based on Liu's definition [50] is "a method that generates a system of algebraic equations for the entire problem domain without consideration of a predefined mesh". This means that the method needs a set of scattered nodes inside the problem domain and on the boundaries of the domain called field nodes. In addition, the relationship between the nodes for the interpolation or approximation of the unknown field variables is not required [49]. Some of these well-known methods are the local point interpolation method (LPIM) [51], the local radial point interpolation method (LRPIM) [52], and the meshless local Petrov-Galerkin method (MLPG) [53,54].
The coupling of LBM and MFree methods has achieved acceptable results in some applications [55][56][57][58]; however, this approach is at an early stage of development and must be improved to address instability in low viscosities and high Re numbers. A key point in the stability and the accuracy of these methods is the collision operator, which is not remarkable in early LB methods such as the BGK model and the multi-relaxation times (MRT) model [59], both of which also violate the principle of Galilean invariance. Therefore, using a suitable collision operator to simulate complex engineering problems accurately has been recognized as a necessity. The advent of the more stable cumulant LBM [60,61] could dramatically improve MFree-LB methods and contribute to their advancement as powerful numerical tools for complex simulations.
The aim of this paper is to study the capability of the local radial point interpolation cumulant LBM (LRPIC-LBM) to simulate the propagation of planar acoustic waves, including the temporal decay of a standing plane wave and the spatial decay of a planar acoustic pulse of Gaussian shape and calculate the deviation from theoretical results, and to determine whether this method might be useful for wind turbine problems. The LB equation is deconstructed into collision and streaming parts. The collision step is performed by means of the cumulant method. The streaming step, which represents a pure advection equation, is discretized first in time using the Lax-Wendroff scheme, and then in space using the local radial point interpolation method (RPIM). Thus, this paper is arranged in the following sections; Section 2 is devoted to a brief summary of our new LB method. Section 2.1 presents the collision part referred to as cumulant collision step. Section 2.2 is devoted to the second most important part of the LBM, the streaming step. The Mfree shape function construction (RPIM shape function), time discretization (Lax-Wendroff scheme), and space discretization (LRPIM) are discussed in Sections 2.2.1-2.2.3. Section 3 reports the results of the LRPIC-LBM. Particularly, the planar standing wave and planar pulse wave propagation results are discussed in Sections 3.1 and 3.2.

The Lattice Boltzmann Method
The lattice Boltzmann method (LBM) is obtained from the kinetic theory of gases. In the LBM, a key point for modeling the momentum distribution in discrete space is the use of probability density functions [22]. The lattice Boltzmann equation without an external force is where f i , Ω i , and c are the particle distribution function, the collision operator, and the lattice speed, respectively. In general, the LBM consists of collision and streaming steps. In the local radial point interpolation cumulant LBM (LRPIC-LBM), the cumulant method is used for the collision part and the local radial point interpolation method (RPIM) is used for the streaming parts.

Collision Step
The cascade method has been proposed [62,63] as a way to overcome the instability problems and modeling artifacts of previous LB methods. However, it is hindered by the effects of lower order moments over higher order moments. The cumulant method, presented by Seeger to solve the Boltzmann Equation [64,65], effectively resolves these issues. Cumulants can be efficiently generated from central moments. The cumulant method for solving the LBM is briefly described in this section.
Using the moment-generating function, M(Ξ, H) = ij f (ξ i , η j )e Ξξ i e Hη j , the moments can be determined without any discontinuity issues as [66] where H, Ξ are the normalized wave numbers. It is a good idea to shift the moment-generating function into the moving reference frame of the fluid to reduce the Galilean invariant issues in the collision process. Therefore, using the central moment-generating function, ı M(Ξ, H) = e −Ξu/c−Hv/c M(Ξ, H), the central moments can be defined as [67] The cumulants are calculated using the logarithm form of the moment-generating function [68] as follows: Each cumulant relaxes with an individual relaxation rate [69].
where c eq ξ m η n are the cumulants of the equilibrium state. The sound speed and the kinematic viscosity are defined as c s = ∆x/( √ 3∆t) and υ = ∆tc 2 s (1/ω − 0.5), respectively.

Streaming Step
To model the streaming step, the second most important part of the LBM, a pure advection equation is normally solved from a Lagrangian approach within uniform structured meshes with CFL numbers equal to one. However, considering the Eulerian perspective for the calculation can effectively resolve this step when meshes are non-uniform and unstructured. The pure advection equation is One alternative is a semi-discrete formulation with time and space derivatives discretized separately. Thus, Equation (7) is discretized using the explicit Lax-Wendroff scheme in time, followed by the local radial point interpolation method (LRPIM) in space. In addition, it is important to approximate the unknown field functions using trial (shape) functions as an approximate solution for the partial differential equation.

MFree Shape Function Construction-Radial Point Interpolation Shape Functions
Radial point interpolation method (RPIM) shape functions were developed to circumvent the singularity problem arising in the point interpolation method (PIM). The RPIM interpolation augmented with polynomials is where R i (x) is a radial basis function (RBF), and p j (x) is monomial in the coordinate space . Parameters n and m are the number of RBFs and polynomial basis functions. Variables a i and b j are time dependent unknown coefficients. It should be noted that the independent variable in RBF R i (x) is the distance between the point of interest x and a node at x i . Different radial basis functions (RBF), and their characteristics have been studied extensively in the meshless RPIM. In this paper, the multi-quadrics (MQ) function is used as where α c = 1.0, q = 1.03, and d c = 3.0.
To determine the n + m unknown coefficients in Equation (8), some specific constraint equations and the Kronecker delta function property are dictated. These constraints are where Thus, the approximation function can be obtained as where F is a vector containing the nodal values of the distribution function and Φ is a vector containing the first n components of theΦ vector where

Semi-Discrete Formulation-Time Discretization
The Taylor series expansion of the particle distributions is where here n refers to the time step. Substituting the time derivatives in terms of the t derivatives up to second order results in the time discretization of Equation (7) based on the Lax-Wendroff scheme,

Semi-Discrete Formulation-Space Discretization
The local radial point interpolation method (LRPIM) was developed to avoid the side effects of using global background cells in the global weak-form. In this method, the numerical integration is performed within the local domain consisting of a set of distributed nodes. The LRPIM is based on the RPIM shape functions with the delta function property. The main advantage of the LRPIM is the excellent interpolation stability of RBFs.
MFree local weak-form methods use the weak form of the problem obtained from the method of weighted residuals (MWR). The weighted residual statement of Equation (16) on the local domain Ω I of point I bounded by Γ I is posed as (17) where w I is the local weight function of node I considered as where r = |x − x i |/d max and d max is the radius of the compact support. Equation (17) is applied to all nodes in the problem domain.
To work with a continuous approximate solution, it is necessary to decrease the differentiation requirements of the unknown in the weighted residual statement by employing integration by parts in Equation (17), where Γ I is the boundary of the local domain Ω I and n β is the unit outward normal vector. Substitution of the approximate solution given in Equations (12) into the weak form given by Equation (19) where M I J and K i,I J are the nodal mass and stiffness matrix, respectively, defined as Thus, the global equation system for all nodes in the entire domain is obtained as where M, K, and f i are the global mass matrix, stiffness matrix, and particle distribution vector, respectively. This system has N equations with N unknowns which is solved separately for each direction.
To numerically evaluate the area and the curve integrals in Equations (21) and (22) the Gauss quadrature scheme is used as follows where n g and n b g are the total number of Gauss points in the quadrature domain and boundaries, w k is the Gauss weight factor for Gauss point x k , J Ω I and J Γ I are the Jacobian matrix for the domain and boundary integrations, respectively.

Results and Discussion
One of the complicated phenomena that has recently received major interest from researchers using the LBM is wind turbine aeroacoustics [32,33,70,71], which can be directly simulated without additional computational cost. In this section, our aim is to demonstrate the standard analysis procedure using the local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM) to predict acoustic properties for benchmark cases. Thus, numerical studies are conducted for the propagation of planar acoustic waves, concentrating on numerical dissipation and dispersion. Considering the total pressure equation as p = p 0 + p , where p 0 is the atmospheric pressure and p is the acoustic pressure, a lossy wave equation is [72] Å 1 + τ s ∂ ∂t where where η is the coefficient of shear viscosity, and η B is the coefficient of bulk viscosity. A quantitative assessment of the method's two setups, including the temporal decay of a standing plane wave in a periodic domain and the spatial decay of a propagating planar acoustic pulse of Gaussian shape for regular and irregular nodal distributions (shown in Figure 1a,b) are discussed below. It should be noted that the default nodal arrangement is a regular distribution; the base units are in the LB system.

Planar Standing Wave
As an initial case study, we performed a temporal analysis on a standing plane acoustic wave in a periodic domain. The dissipation and dispersion relations based on the temporal analysis of Equation (26) are [25] where k is the wave number. The assumptions for this set-up are presented in Table 1. They were chosen as in reference [26] for ease of comparison. Variables The acoustic pressure at time t is [26] p (x, It should be noted that the results of the temporal analysis can be considered as a function of the number of points per wavelength N ppw = λ/∆x or the non-dimensional wave-number k∆x = 2π/N ppw . In accordance with the concepts discussed in [26] to study the accuracy of the propagation of waves using the local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM), various viscosities and different resolutions (i.e., the number of points per wavelength) were studied. Figure 2 shows the acoustic pressure time history for ν = 1.0 × 10 −2 ( ∆x 2 ∆t ) with N ppw = 12 points per wavelength. The analytical result is represented by a solid black line, whereas the cumulant LBM and the LRPIC-LBM are shown with a red dashed line and a blue five-pointed star, respectively. The deviations of the numerical phase speed and the temporal dissipation rate from the theoretical values are less than one percent for the BGK LBM [25] and the cumulant LBM [26] for resolutions lower than 12 points per wavelength. However, the LRPIC-LBM exhibits even better behavior in predicting the acoustic pressure of the analytical values. The acoustic pressure time history for ν = 1.0 × 10 −4 ( ∆x 2 ∆t ) with N ppw = 12 points per wavelength is presented in Figure 3. It shows the comparison between the analytical solution, the cumulant LBM and the LRPIC-LBM numerical solution at low viscosities. One of the most problematic issues in the BGK LBM is the low viscosity limit, which makes the solution unstable. However, the cumulant LBM, with phase speed and temporal dissipation rate errors of less than 1 percent, as in reference [26], did not present difficulties at low viscosity. Although both approaches were successful, the new method exhibited better performance, closely following the theoretical result for the same viscosity value. An important characteristic highlighted by some researchers [25,26] is that these numerical deviations are only a function of N ppw and are independent of other parameters such as frequency and viscosity. They found that the errors of the BGK [25] and the cumulant LBM [26] are about 7 percent for N ppw = 4. However, the acoustic pressure time history for ν = 1.0 × 10 −2 ( ∆x 2 ∆t ) with N ppw = 4 points per wavelength illustrated in Figure 4 for the analytical solution, the cumulant LBM and the LRPIC-LBM reveals that the deviation is less than 2 percent for the LRPIC-LBM, with ∆t = 0.25. Thus, the LRPIC-LBM is much more successful in predicting theoretical results even with a low number of points per wavelength. The choice of the time step size has an influence on the accuracy and stability of the solution. A smaller time step leads to more accurate results, especially in a hyperbolic partial differential Equation (PDE). To minimize the phase speed and dissipation rate errors, the time step was reduced. Figure 5 shows the acoustic pressure time history for ν = 1.0 × 10 −2 ( ∆x 2 ∆t ) with N ppw = 4 points per wavelength and ∆t = 0.1. The LRPI-CLBM replicates the analytical results with negligible errors. Therefore, this method makes it possible to predict wave motion more accurately with no dependency on the number of points per wavelength N ppw . Although the results of the propagation of acoustic waves for regular nodal distributions were good, it is important to adequately estimate accuracy when considering irregular nodes (Figure 1b). The acoustic pressure time history for ν = 1.0 × 10 −2 ( ∆x 2 ∆t ) with N ppw = 12 points per wavelength is presented in Figure 6. The results show that the local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM) with irregular nodal distributions again very closely reproduced the analytical acoustic pressure.

Planar Pulse Wave
Next, we studied a planar pulse wave by replacing the plane wave with a Gaussian shape planar pulse, initially located at the center of the domain. The dissipation and dispersion relationships derived from the spatial analysis of Equation (26) are [25] The variables and parameters for this case are presented in Table 2. A planar pulse emits from the origin throughout the domain, where periodic boundary conditions are imposed.
To assess the accuracy of the LRPIC-LBM in simulating the pulse wave emission, we proceeded as in reference [26]. Different viscosities and resolutions (which are related to the choice of σ) were studied. Figure 7 depicts the acoustic pressure time history for ν = 1.0 × 10 −2 ( ∆x 2 ∆t ) with σ = 0.1. The cumulant LBM results (dashed red line) are compared to the LRPIC-LBM (solid black line). As stated in [25,26] the intensity loss at any location is proportional to the distance of propagation, regardless of the precise location. Thus, the data were extracted from the center of the domain, and at 5, 11, and 17 nodes apart from the center. For resolutions up to σ = 0.1, deviations between the cumulant LBM and the LRPIC-LBM results were less than 1 percent. The acoustic pressure time history for ν = 1.0 × 10 −4 ( ∆x 2 ∆t ) with σ = 0.1 is shown in Figure 8. It presents the comparison between the cumulant LBM and the LRPIC-LBM at low viscosity. While the standard BGK LBM creates instabilities and noisy results, less than 1 percent deviation was found between the cumulant LBM and this method, with stable behavior at low viscosities. In addition, as in the case of the temporal results shown in Figure 3, reducing the viscosity in the considered range did not substantially impact the results.
Like the parameter N ppw introduced in the temporal analysis of the first set up, σ is another effective parameter used in spatial analysis which affects the accuracy of the LBM results.The acoustic pressure time history is depicted in Figure 9 for ν = 1.0 × 10 −2 ( ∆x 2 ∆t ), with σ = 0.06. It shows that the deviations between the cumulant LBM and the LRPIC-LBM results increase after the reduction of σ. The wiggling observed in the results of the cumulant LBM is due to the reduction of the number of nodes inside the pulse, which brings the accuracy of the results into question. However, the LRPIC-LBM graphs are smooth and unaffected by the lesser number of nodes. To better assess the accuracy of the LRPIC-LBM compared to the cumulant LBM, the data shown in Figure 9 were compared with the analytical results. The Fourier transform of the pressure time history of the passing wave yields the Fourier coefficient of pressure which is the solution to Equation (26) [25]. Figure 10 shows the ratio " p (x, ω)/ " p (x, 0) as a function of angular frequency for the analytical solution. This figure shows that the LRPIC-LBM is more successful than the cumulant LBM at predicting theoretical results with σ = 0.06. Thus, this method can model wave propagation more accurately even at smaller resolutions.

Conclusions
This paper presents a numerical study of the propagation of acoustic waves, one of the challenging issues occurring in wind turbine simulations, that includes the temporal decay of a standing plane wave and the spatial decay of a planar acoustic pulse, using the local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM). The LB equation is divided into collision, modeled by the cumulant method, and streaming, discretized in time and space, using the Lax-Wendroff scheme and the local radial point interpolation method (RPIM). The LRPIC-LBM results were compared with the cumulant LBM and the analytical solutions. Both methods showed a similar acoustic pressure time history, and the deviations for the phase speed and the dissipation rate were minor for high number of points per wavelength N ppw and resolution σ. In addition, they showed that reduced viscosity does not affect the stability of either LB method due to the intrinsic characteristics of the cumulant method. Unlike LB methods such as the BGK LBM and the cumulant LBM, the time history for the acoustic pressure and the phase speed and dissipation rate predicted by LRPIC-LBM showed considerably smaller errors for low N ppw and σ due to the construction of the meshless method itself. Moreover, the LRPIC-LBM with irregular nodal distributions reproduces the same propagation of acoustic waves obtained with regular nodal distributions.
In summary, the freedom to scatter nodes based on problem conditions and the occurrence of sharp gradients plus the accuracy obtained from the RPIM along with the good stability and simplicity achieved by the cumulant LBM may provide an adequate platform with which to model wind turbine problems.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: