Aeroacoustic source term computation based on radial basis functions

SUMMARY In low Mach number aeroacoustics, the known disparity of length scales makes it possible to apply well‐suited simulation models using different meshes for flow and acoustics. The workflow of these hybrid methodologies include performing an unsteady flow simulation, computing the acoustic sources, and simulating the acoustic field. Therefore, hybrid methods seek for robust and flexible procedures, providing a conservative mesh to mesh interpolation of the sources while ensuring high computational efficiency. We propose a highly specialized radial basis function interpolation for the challenges during hybrid simulations. First, the computationally efficient local radial basis function interpolation in conjunction with a connectivity‐based neighbor search technique is presented. Second, we discuss the computation of spatial derivatives based on radial basis functions. These derivatives are computed in a local‐global approach, using a Gaussian kernel on local point stencils. Third, radial basis function interpolation and derivatives are used to compute complex aeroacoustic source terms. These ingredients are necessary to provide flexible source term calculations that robustly connect flow and acoustics. Finally, the capabilities of the presented approach are shown in a numerical experiment with a co‐rotating vortex pair.

geodesy, mapping, signal processing, digital terrain modeling, and hydrology. 2 In the field of aeroelasticity, fluid-structure interaction algorithms based on mesh-less radial basis function (RBF) interpolation are proposed as a computational efficient coupling strategy. [3][4][5] The coupling of nonconforming grids with focus on computational performance is of great importance in aeroelasticity. 6 Today, efficient interpolation techniques are especially important when dealing with computational aeroacoustics. Since the beginning of CAA, hybrid methods have been established as the most practical for fast and accurate aeroacoustic computations at low Mach numbers. The workflow of hybrid aeroacoustics involves three steps (see Figure 1): (i) perform unsteady flow computations on a restricted subdomain; (ii) compute the acoustic sources; (iii) simulate the acoustic field. For low Mach number applications, a key challenge in computational aeroacoustics (CAA) is the huge disparity of scales between flow structures and audible acoustic wavelengths. The scaling between the acoustic wavelength and the characteristic length l of a vortex is given by where Ma is the Mach number. Utilizing this disparity of scales, aeroacoustic analogies or perturbation techniques separate the flow computation from the acoustic computation (eg, References 7,8). Therefore, the sizes of the two computational grids are in general quite different. For both physical fields, the individual optimal computational grid achieves the highest accuracy and the two grids differ according to the modeling criteria. On the one hand, the flow grid resolves boundary layers and is mostly coarsened toward outflow boundaries to dissipate vortices. On the other hand, the acoustic grid transports waves and therefore needs a uniform grid size all over the computational domain. In order to couple the different meshes, an interpolation scheme is necessary that satisfies the fundamental requirement of hybrid CAA: An accurate data transfer from the flow to the acoustic grid to minimize interpolation errors while conserving the energy. To cope with this task, different interpolation strategies can be applied, starting from low complexity nearest neighbor interpolation to complex source-term computation models with volume intersections between flow and acoustic grid. The simple nearest neighbor interpolation fails to compute the acoustic sources accurately (eg, References 9,10). In Reference 9, the aeroacoustic source term is computed by summing the contributions of all flow cells belonging to finite elements that surround an acoustic finite element node. It is assumed that the flow quantities used for the acoustic source-term computation are constant over each flow cell. This approach has also been used in Reference 11 for three-dimensional problems, where a grid dependency has been observed resulting in a too low sound pressure level over the whole frequency spectrum. A fully conservative approach has been used in Reference 12, where the acoustic sources within the finite element formulation are first computed on the fine flow grid. These so-called nodal loads are then interpolated by a conservative scheme to the acoustic grid. This approach is accurate in cases where the flow grid is much finer than the acoustic grid but fails in cases where the flow grid is coarser than the acoustic grid. As a solution to this problem, we derived a cut cell approach and successfully applied it to the aeroacoustic computation of an axial fan. 13 Similar investigations have been performed in Reference 10, where for both the flow and the acoustic field a finite volume scheme has been used. However, most of these methods fail if the computational fluid dynamics (CFD) simulation only provides primary variables, such as the pressure p or the velocity u, since the aeroacoustic sources are often nontrivial combinations of these primary variables (see Figure 1).

F I G U R E 1
The hybrid aeroacoustic workflow consists of three main computational parts, where the presented radial basis function source-term computation is a crucial part if the CFD solver provides only primary results, such as the pressure p or the velocity u. Here, a prime denotes the fluctuating part of a physical quantity and a bar over a physical quantity its time average [Colour figure can be viewed at wileyonlinelibrary.com] Desirably, an accurate, flexible, and conservative coupling scheme ensures a rigorous connection between fluid dynamics and acoustics within a hybrid aeroacoustic simulation. The properties of the desired coupling scheme, or interpolation, are summed up as follows: • An accurate and fast interpolation technique.
• An interpolation technique that handles special grids, for example, grids to resolve boundary layers.
• A method to compute accurate derivatives of the primary flow variables, such as pressure p, velocity u, density , temperature T, and entropy s.
• A flexible algorithm that can be integrated into a standard product development cycle.
• A flexible algorithm that assembles different hybrid aeroacoustic source terms, for example, the divergence of the Lamb vector ∇ ⋅ (((∇ × u) × u) ′ ) or the divergence of an entropy source ∇ ⋅ (T ′ ∇s − s ′ ∇T). 14 • A conservative algorithm that transfers the desired amount of energy, defined by the aeroacoustic sources, from the flow discretization to the mesh of the acoustic simulation (see eg, Reference 13;) is not discussed here.
In this paper, we propose to compute the aeroacoustic sources directly from the primary CFD variables by applying an interpolation scheme based on RBFs in conjunction with RBF derivatives. We show the ideal setup for the algorithms, the choice of the kernel function and propose a selection of discretization natural neighbors for the computation based on the connectivity. This patch search technique guarantees the resolution of typical flow structures. The application of local RBFs provides promising capabilities in terms of computational efficiency, known from nearest neighbor algorithms. Furthermore, the computation of RBF derivatives can be carried out elegantly and accurately with a local-global approach.
The rest of the paper is organized as follows: In Section 2, we provide a general overview about the selected interpolation scheme, discuss convergence based on the parameters, the grid characteristics, the neighbor search, and typical fluid dynamic effects. Boundary layer interpolation motivates the proposed neighbor search technique and the superior performance is shown in a convergence study. Furthermore, the integration of constraints is discussed. Next, in Section 3, we discuss our approach to compute spatial derivatives based on local-global RBFs, where we investigate the sensitivity of the scaling parameter based on verification benchmarks. The derivatives are judged by analytic representations of flow structures, such as a vortical distribution, a sinusoidal shock, and a boundary layer profile. Finally, in Section 5 we conclude the findings and give additional remarks for future research.

INTERPOLATION SCHEMES
RBFs are scalar-valued multivariant functions R d  → R that depend on the radial measure r in terms of the Euclidean norm r = ||x − z|| 2 between a scattered data (source) point x and the evaluation (target) point z (see Figure 2). Over the years, application-driven development of RBFs has modified function type, shape, and relevant support. The most prominent RBF kernels are the Gaussian, the multiquadric, and the inverse multiquadric. 15 At this point, a clear distinction between local and global methods should be drawn. Global methods mainly using the Gaussian kernel are well suited for theoretical analyses since it can be shown that the solution converges to the exact solution under certain conditions. 16 But on the downside, the condition number of the global interpolation matrix tends to infinity. There are some approaches to circumvent this issue 17 using a QR-decomposition or a Taylor-expansion but they are computationally more expensive and less parallelizable than local ones.
With increasing number of ill-placed datapoints in datasets, the demand for RBFs with compact local support 18 increases. Local representations reduce computational time when interpolating data and increase the condition number of each local subproblem. Instead of inverting one ill-conditioned global system matrix, many small well-conditioned stencil-matrices are inverted. Using Wendland's compactly supported RBFs together with a modified Shepard's method as presented in Reference 19, high computational efficiency for large system can be achieved.
An important aspect of this work is the application of RBF interpolation schemes to scattered data distributions with bad quality, such as extremely anisotropic point distribution, for example, from boundary layer flow. We use a local collocation method that represents the data exactly in the prescribed point. In this sense no algorithmic uncertainty is present at the individual, deterministic, scattered data point. To handle anisotropic point distributions, a special search procedure is used to find the optimal local point distribution for interpolation, which will be presented in Section 2.2.

RBF interpolation
As mentioned above, the main purpose of our method is to interpolate flow results, such as velocity u or vorticity , from a CFD mesh to an acoustic finite element mesh (FEM). Since hybrid aeroacoustics deals with a large number of unknowns and unfortunate data distributions, for example, in boundary layers, we focus on a parallelizable algorithm to achieve high computational efficiency. In this work, the local Wendland kernel Φ (1) together with a modified Shepard's method was chosen, similar to Reference 19. This is the most promising approach since it is both fast and capable of handling boundary layer meshes if the chosen kernel is combined with the patch search technique presented in Section 2.2.
The scaled Wendland kernel 19 reads where x ∈ R d is the location of a scattered data point. For all scattered data points the associated indices s ∈ 1, 2, … , N s are defined, with N s as the number of scattered data points, also called source points. The point z ∈ R d is the point at which the interpolation is evaluated. These points are called target points and form a set of target point indices t ∈ 1, 2, … , N t , with N t as the number of target points. The parameter is responsible for the scaling of the compact support (1), which is especially important for bad source point distributions resulting in a high condition number of the interpolation matrices. For our investigations, the parameter is chosen to be with I q as the indices of the N q neighbors of the target point z and define the patch set. It should be noted that with this RBF interpolation, a higher order of accuracy can be achieved by "flattening" the compact support of the kernel (see Section 3.1) or increasing the order of the local interpolation kernel.

Local interpolation method
Following the approach of Lazzaro, 19 two different scattered data patches around the target point at which the interpolant is evaluated are introduced. The set X q = {x i ∈ x, i ∈ I q }, with I q as the indices of the N q neighbors of the target point z, defines the patch set. The set X w = {x i ∈ x, i ∈ I w }, with I w as the indices of the N w influence points, defines the influence of the different patches. The influence radius r W k = max k∈I w {r k }, with r k = ||x k − z|| 2 , is defined to be the maximum distance of a target point in X w and N w ⊂ N q . The interpolant is given by where R k (z) defines the local interpolation system and W k the weight function. The local interpolation system is given by the algebraic system which has to be initially solved for the N q temporary (unscaled) interpolation weights c j , where R k (x k ) denotes the scattered data (field that is interpolated) at the patch source point x k . The weight function W k of the modified interpolant s(x) is chosen to be with the exponent p as a measure of locality. The bigger p, the more local the approach making it less accurate but capable of resolving stronger gradients. Since the weights W k (x) constitute a partition of unity, the following holds: Then the interpolant s(x) can be evaluated at every target point z. The algorithm of the local RBF interpolation is written in pseudo code (Algorithm 1).

Algorithm 1. Local RBF Interpolation
for k in t do 4: //I q holds indices of N q source points in patch 5: x k = GetCoordsOfSrcPointsForTrgPoint(k) 6: //get scattered data at patch source point 7: Solve the N q × N q -system for the patch coefficients c k , Eq. (4) 8: i w = I q [influence-point] //get index of influence point 9: r W k = GetDistOfPoint[i w ] //distance between target point and i w -th source point 10: Compute weights for every source point in patch k: W k , Eq. (5) 11: Compute interpolant s(z k ) at target point z k , Eq. (3) 12:

end for 13: end function
Proof of convergence of the interpolation. As described by Lazzaro in Reference 19, an upper bound for the global approximation error E(x) reads with e k (x) = ||f (x) − R k (x)|| as the approximation error of the interpolant R k (x) relative to the local interpolation system. Additionally, Madych has proven in Reference 16 that the following holds if the distance h between the source points is sufficiently small: where C 1 , C 2 ∈ R + are constants. Thus for the shape parameter of the Wendland kernel (1) → ∞, the approximation error of the local interpolation system e k (x) goes to zero and hence from Equation (7) it can be concluded that the same is true for the global approximation error E(x). ▪

On the optimal choice of neighbors
An important part of the interpolation is the choice of source points x i in the patch X q to avoid unphysical artifacts in the interpolated field. This is especially important in boundary layers, where cells or elements are distorted or stretched with high aspect ratios. A common approach to obtain a specific amount of nearest neighbors is to use a kd-tree search, for example, from computational geometry algorithms library 20 or fast library for approximate nearest neighbors. 21 These fast kd-tree algorithms have to be provided with the coordinate of the target point for which the patch is defined and all coordinates of the source points. By doing this, the underlying mesh structure (connectivity) is ignored and the data is represented on discrete scattered points (see Figure 3A). Thus, important information to improve the quality of the interpolation might be missed. To understand the shortcoming of this approach, a simple example of a boundary layer flow in a channel is provided, where the cell data from a CFD simulation is interpolated to nodes of a different FEM mesh. The interpolation from a CFD Finite Volume mesh to an FEM mesh is our standard procedure; however, the findings are also true for other discretization types without loss of generality.
The interesting part is the area close to the wall, where a mesh refinement is necessary to properly resolve the boundary layer. The resulting patch of elements emerging from a kd-tree search for two example nodes on the boundary layer is shown in Figure 3A. There we can see the unfavorable choice of source points (black) for different target points (red circle). As shown in Figure 3B, this patch choice results in a poor interpolation in streamwise direction, which leads to numerical artifacts. The results become even worse if we provide a heavily distorted mesh. There we can observe that the patches only contain vertical neighbors and thus the patches have no connection in streamwise direction. A common approach to circumvent this issue is using anisotropic search algorithms and anisotropic kernel functions. In contrast to the herein described connectivity-based neighbor search, additional tuning parameters have to be determined for local mesh features.
However, if we use a connectivity preserving mesh-based neighbor search, the patch looks as displayed in Figure 4A, which results in an even spatial distribution and good interpolation results (see Figure 4B). This modified patch algorithm searches for the next directly connected neighbor elements over the common global nodes. By applying a "layer strategy," the neighbors of the direct neighbors are also taken into account. This search can be carried out recursively in a layered manner until the preferred number of elements is reached. In the illustrated example, only a one-layer level search was applied.
Another drawback of a coordinate based kd-tree search is the interpolation along slender wing profiles and structures. It is possible that a kd-tree search at the upper side of the wing includes points from the lower side in the patch, which would be completely incorrect. In this case, the modified patch search based on the connectivity performs well since elements on the upper side have no direct connection to the elements on the lower side of the slender airfoil.
With the given interpolation formulation it is possible to impose exact boundary conditions, for example, a no-penetration condition u ⋅ n = 0, where the corresponding entries of the interpolant can simply be set equal to zero.

Convergence of interpolation
To analyze the convergence of the RBF interpolation and to quantify the more or less heuristic observations of correct choice of patches, we prescribe analytic functions on a unit cube Ω ∈ [0, 1] 3 with different node distributions and observe the L 2 -error after interpolating to a different mesh with equidistant node distribution.
The first analytic function represents the transition of different vortex structures in shear layers (see Figure 5A). This smooth function is used to investigate the convergence and accuracy of the RBF interpolation since the majority of pressure fields encountered in low Mach number aeroacoustics are smooth. The second analytic function is a representative of a continuous sinusoidal shock, which is similar to the numerical solution of a shock using a dissipative flow solver 22 (see Figure 6A). This function is used to test the capabilities of RBF interpolation to handle strong gradients and sharp transitions. These analytic functions were prescribed on a unit cube Ω ∈ [0, 1] 3 consisting of equidistant hexahedrons, and tetrahedron elements. Then the analytic functions were interpolated to a different mesh consisting of equidistant hexahedrons with discretization distance Δz using a nearest neighbor search or the modified patch search, respectively. Figures 5B and 6B show the convergence of the RBF interpolation with respect to the source discretization distance Δx, being the edge length of an element, for the first and second analytic function. Obviously, there is an influence resulting from the source data density, and the RBF interpolation performs well for both the nearest neighbor search and the modified patch search on all of the different mesh types. The larger L 2 -error of the tetrahedron elements is due to an unsymmetric interpolation point distribution around the target point, which is suboptimal for RBF and leads to higher errors. However, the order of convergence remains the same.
In order to further analyze the capability of the RBF interpolation to handle boundary layer data in combination with the different choices of neighbors, we use the function ) .
The analytic function f 3 (x) roughly describes the flow velocity component in inflow direction in a boundary layer (see Figure 7A). Here, the wall that is the reason for the boundary layer to develop is located at y = 0. As mentioned above, a CFD mesh to properly resolve the boundary layer is characterized by cells with a big aspect ratio in the region adjacent to the wall where the largest gradients occur. For our investigations, we used a mesh consisting of regular hexahedrons with different aspect ratios in the range of 2 to 512 for the cells at y = 0 (typical CFD boundary layer mesh 23 ). With growing distance y the mesh gradually becomes coarser. In order to realize different aspect ratios in the first cell row, the cell height is varied while the number of nodes on the edges stays the same.
As shown in Figure 7B, the kd-tree-based nearest neighbor search interpolation causes increasing errors with increasing aspect ratio. Hence, the interpolation error for large aspect ratios limits the application of the nearest neighbor search. However, the previously described patch search leads to better results at large aspect ratios where it significantly outperforms the kd-tree-based nearest neighbor search. Ideally, the L 2 -error is independent of the aspect ratio. One can see that this requirement is approximately met using patch search.

DERIVATIVES BASED ON RBFS
Although the procedure to compute derivatives based on RBFs is similar to the RBF interpolation, there is one important difference. For derivatives, there is only one source patch set X q = {x i ∈ x, i ∈ I q }, with I q as the indices of the N q neighbors of the target point z. To explain the derivative of a given function s(z) at target point z, let us start using matrix notation. We first need to interpolate the values at discrete datapoints x to determine the interpolation weights . Thus we need to form the interpolant s(x) and solving this system for the weights leads to Since we want to apply the differential operator [ ] to the interpolant, we may write and inserting the weights from Equation (13) yields Now we can use the fact that = −1 (x) ⋅ s does not depend on z and we can pull it out of the operator's range and we can rewrite more compactly in the following way: The process from above can now be transformed to a more implementation-friendly form, based on Reference 24, where we first have to apply the spatial derivative operator  to the RBF and evaluate it, as presented in Equation (18), where x k and x l are the coordinates of two source points with indices k, l ∈ I q .
whereas the derivative of the Gaussian kernel was computed analytically. This system has to be inverted in order to obtain the weights c l for the derivatives at the source points l ∈ I q . Then the derivative of the given function s(z) can be evaluated by where R l = s(x l ) ∈ R N q ×1 is the scattered data value of the given function at the source point with index l. Algorithm 2 shows how to calculate derivatives of scattered data distributions. Vector-valued datasets can be interpreted as a tuple of scalar components ) . (20) Then R k is a matrix R N q ×d with d as the spatial dimension and evaluating the last expression of Algorithm 2 leads to a Jacobian matrix (21), where each row i contains the derivatives of the scattered data vector.  x k = GetCoordsOfSrcPointsOfTrgPoint(t) 6: for k in I q do 7: //Compute the weight for source point k, Eq. 18 10: end for 11: R t = GetSourceData(I q ) //Get the scattered data at patch source point 12: s(z t ) = c ⋅ R t //Derivatives s(z t ) at z t , Eq. 19 13: end for 14: end function To show how to form the desired differential operator from the spatial differentiation matrix and to outline the workflow, a simple example is presented. Assume we have a scattered data vector field (s 1 , s 2 , s 3 ) T (z) ∈ R 3 . Then the rows of L i Φ contain the spatial derivatives of the three-dimensional RBF in the three spatial directions . Additionally, the derivative coefficients c[i] for source point i are three-dimensional vectors. The final evaluation S(z) = c ⋅ R(z) leads to the following 3 × 3 matrix

Kernel scaling
As mentioned above, we use a Gaussian kernel in a local approach. The local kernel has no numerical stability problem in its basic form and works well on regular grids. But as soon as the mesh becomes distorted (eg, boundary layer) the results worsen. An improvement can be achieved by scaling the Gaussian ansatz function in the following way: We found that the Gaussian kernel theoretically leads to high accuracy (in theory we can differentiate up to any given order) but the stability problem, similar to the global interpolation method, restricts the maximum flatness . Roberts proposed a method using a Taylor's expansion in the flatness parameter to circumvent this issue. 25 However, there are optimization strategies 15 to find optimal values of the flatness parameter .
It should be emphasized here that the parameter is computed differently depending on whether an interpolation or derivation procedure is carried out. Our idea to find an optimal shape parameter for the RBF derivative procedure is to derive a relationship (a confident initial guess) for the parameter with respect to the value r 12 as the smallest distance between two points in the stencil of the current target point. The basic approach is to minimize the following objective functional with the L 2 -error E, the Lagrange multiplier , the condition number , and the limit of the condition number limit where the system becomes sensitive to floating-point round-off errors. The L 2 -error represents the deviation of the basis function from being one (absolutely flat, → ∞) on the whole domain. This error measure can simply be obtained by integrating The challenging task is to find the influence of the shape parameter on the condition number of the resulting linear system. In general, the condition number of a matrix is defined as the ratio of the largest to the smallest eigenvalue: Therefore, the smallest eigenvalue must be nonzero in order to stay invertible. The minimization is trivial since the functional decreases monotonically in while a critical condition number limit guarantees the invertibility of the matrix. With this constraint, our confident initial guess initial is initial ∼ 1∕r 12 (27) and from this starting point the shape parameter is iterated to the limit of the matrix invertibility, depending on the inversion algorithm. As proven multiple times in the literature dealing with Gaussian RBF kernels (eg, in Reference 19), the L 2 -error decreases with increasing shape parameter . However, this is true only up to a certain value, where numerical instabilities tend to increase the error again.

Additional notes
During the computation, there is an easy way to assess how accurate the derivative will be. Looking at the last expression of Algorithm 2, where we multiply the weight matrix c ∈ R d×N q and the source data matrix R ∈ R N q ×d , we can state that the rowsum of each row must be zero. This can be seen easily by considering a constant field, where all entries are equal R ij = const. and the resulting derivative must be zero (derivative of a constant function) S i (z) = ∑N q j=0 c ij R j = 0. For nontrivial coefficients, this only holds if the rowsum is zero ∑N q j=0 c ij = 0 for i = 1, … , d. This rowsum condition is used for optimal neighbor consideration.
During the investigations, the local Wendland kernel Φ 2,0 = (1 − ||x − z|| 2 ) 2 works perfectly for generic test cases with moderately distorted elements. But as soon as the Wendland kernel is applied to more sophisticated problems with boundary layers, severe numerical artifacts are observed, especially when computing the divergence. Therefore a local approach with a Gaussian kernel (22) is used, which produces accurate results for a wide range of mesh discrepancies.

Convergence of derivatives
To analyze the convergence and accuracy of the RBF derivatives, we again prescribe our analytic functions Equations (9) to (11) on a unit cube Ω ∈ [0, 1] 3 with different node distributions using the patch search algorithm and observe the L 2 -error after computing the gradient of the field. Figure 8A shows the convergence of the RBF derivatives with respect to the characteristic discretization distance Δx for the first and second analytic function. Similar to the RBF interpolation, there is an influence resulting from the data density. Although the RBF derivatives are accurate for both analytic functions, it is obvious that larger gradients (as encountered in the second analytic function) lead to larger deviations.
Additionally, the convergence of the derivatives in boundary layers was analyzed considering the analytic function (11). This analytic function was again prescribed on a mesh consisting of regular hexahedrons with different aspect ratios in the range of 2 to 128 for the cells in the first row.
As shown in Figure 8B, the L 2 -error decreases with increasing number of nodes in the region where the largest gradients occur (next to the wall). However, since the number of nodes on the edges perpendicular to the wall stays the same while the height of the cells is varied, a point is reached eventually where the nodes are not distributed well enough over the entire region of the largest gradients. This is the reason for the flattening of the curves for hexahedrons with large aspect ratios.

NUMERICAL EXAMPLE
The rotating vortex pair has been frequently used to determine the capabilities of aeroacoustic methodologies. [26][27][28][29] This arrangement has the nature of a quadrupolar sound field. Figure 9 illustrates the configuration of the vortex pairs. Both vortices are delta distributions and oppose each other at a distance of 2r 0 . The strength of each vortex is characterized by the circulation intensity Γ. The vortices rotate around the origin with a period of T = 8 2 r 2 0 ∕Γ imposing an angular rotating speed r = Γ∕(4 r 2 0 ) ⋅ e 3 = r ⋅ e 3 . Each vortex convects the other vortex by a velocity of u = Γ∕(4 r 0 ) ⋅ e t , where e t is the unit vector in tangential direction. The Mach number in the circumferential direction is given by M = u ∕c = Γ∕(4 r 0 c). As already mentioned, the hydrodynamic field of the spinning vortex pair is expressed by a rotating quadrupole field. The potential flow theory can be used to determine the fundamental solution of the spinning vortex pair in terms of the complex flow potential function. In doing so, we introduce the transformation from Cartesian coordinates (x 1 , x 2 ) to the complex plane with the complex coordinate z = rexp i = x 1 + ix 2 . The location of each vortex over time t is defined by b = r 0 exp i t . Using the definitions, we express the incompressible, inviscid flow potential (z, t) as From this field we can derive the primary variables representing the fluid dynamic quantities. Based on the pressure and the velocity, aeroacoustic sources are derived. The incompressible velocity field u = (u 1 , u 2 ) T of the spinning vortex pair is obtained by differentiating Equation (28) with respect to the complex coordinate z The fluid dynamic pressure p ic is obtained by applying the unsteady form of Bernoulli's principle Müller and Obermeier 30 derived an analytic solution of the acoustic far-field, based on matched asymptotic expansion of the potential solution. Starting from the solution in form of a complex potential (28), matching the inner and the outer solution yields the far-field of the acoustic pressure fluctuation p ′ of the co-rotating vortex pair In Equation (31) k = ∕c denotes the wave number, J 2 (⋆) the second-order Bessel function of first kind and Y 2 (⋆) the second kind. This analytic pressure fluctuation p ′ is used to verify both aeroacoustic methodologies, that is, using the divergence of Lamb vector and the first substantial time derivative of the incompressible pressure as acoustic source terms. It should be emphasized that this fluctuating pressure p ′ is not equal to the acoustic pressure p a ; however, p ′ → p a holds in the far-field.

Numerical investigation
In order to judge the computational procedure, [26][27][28][29] we compare the analytic solution of the test case with a numerically obtained solution based on the simulation procedure that is illustrated in Figure 1. This analytic validation has the benefit of mitigating experimentally induced errors, reducing the overall modeling uncertainty of the aeroacoustic model, and simply avoiding errors of a numerical benchmark routine. We simulated the co-rotating vortex pair on a stationary grid with moving sources induced by the vortical structures. An unstructured mesh is used to discretize the computational domain. In the source region, a characteristic mesh size of h ≈ 90 cm is used. Each vortex distribution Γ (z − b) is approximated by a continuous multivariant normal distribution with equivalent circulation Γ and an isotropic variance of 2 = 0.05 m 2 .
We assume that the analytic field (u, p ic ) is represented on the flow grid and is based on a circulation strength of Γ = 2 m 2 ∕s and a distance of 2r 0 = 2 m between the vortices. The angular rotation induced by the vortices is r = 0.5∕s, the speed of sound c = √ 10 m/s and density 0 = 1 kg/m 3 . Using this flow field, we apply the source term computation procedure and simulate the sound which is compared to the analytic solution in the far-field (31). We verify the interpolation methodology by two hybrid aeroacoustic simulations of the co-rotating vortex pair. The first simulation is based on vortex sound and the divergence of the Lamb vector as source term 31 The second verification example is based on the perturbed convective wave equation (PCWE) with the substantial derivative of the incompressible pressure as aeroacoustic excitation 13 1 In Equation (33), a denotes the acoustic velocity potential and D Dt = t + u ⋅ ∇ the substantial derivative. The acoustic pressure can be computed as p a = 0 D a Dt . Based on the computational methodology of RBFs, we propose the following simulation workflow for the computation of the aeroacoustic source terms. First, the primary fields incompressible pressure p ic , flow velocity u ic , and vorticity are computed in a flow simulation on the CFD mesh. Second, if required, derivatives are computed by the RBF framework and the source term is assembled on the CFD mesh. Naturally, the RBF framework is carried out on the discretization of the flow simulation, since this discretization already resolves characteristic flow structures. As a third step, a conservative integration 13 of the source terms ensures a nondissipative calculation of the aeroacoustic finite element sources that are applied to the finite element simulation. Finally, the acoustic simulation is carried out, typically on a much coarser CAA mesh.
The source term of Equation (32) involves the cross-product to determine the Lamb vector × u and employing RBF derivatives yields the divergence of the Lamb vector ∇ ⋅ ( × u). Of course, the presented framework is capable of calculating the vorticity as the curl of the velocity field = ∇ × u, with the source term ∇ ⋅ ((∇ × u) × u).
For the PCWE (33), the proposed workflow is as follows. The substantial derivative of the incompressible pressure is computed from primary CFD quantities. With the use of RBF derivatives, the gradient of the incompressible pressure is evaluated and the inner product with the velocity field is computed. We utilize a fifth-order backward differencing scheme for the time derivative of the pressure. Then the source term is assembled for the acoustic simulation.
For both acoustic simulations, a Newmark scheme with a time step size of Δt = 0.09 s is applied for the time discretization. The computational domain (260x260m) is structured in three separate conforming subdomains. A source domain consisting of a disk with a diameter of 5r 0 (relatively high resolved unstructured triangular mesh, h ≈ 90 cm) is embedded in a propagation region, in which we gradually increase the unstructured mesh size. The propagation region is surrounded by a structured perfectly matched layer region, 32 which absorbs the radiating waves. The wave length = 2 k ≈ 20 m is resolved (in the propagation region) with approximately 20 linear finite elements per wavelength.  Figure 10A shows the acoustic field of co-rotating vortex pair, with the Lamb vector as source term. The acoustic field with the substantial derivative as source term is depicted in Figure 10B. In the acoustic near-field, differences occur due to the different solution quantities in the two aeroacoustic formulations (fluctuating pressure p ′ with the divergence of the Lamb vector as source term and acoustic pressure p a for the PCWE). Both results have the characteristic radiation pattern of the co-rotating vortex pair.
The steepest descent of the Gaussian distribution is discretized by five linear triangular elements over 2 . This coarse approximation of the vortical distribution shows the robustness of the RBF derivatives. A mesh refinement of the source region (see Figure 11) causes no significant increase in accuracy compared to the analytic solution. As depicted in Figure 11, one can clearly see the good accordance of numerical and analytic solution, even for the coarse grid.

SUMMARY AND CONCLUSION
In this work, we have presented how to use RBFs in a local formulation for the interpolation of scattered data using a Wendland kernel in combination with a modified Shepard's method. Since hybrid aeroacoustics deals with a large number of unknowns and unfortunate data distributions, for example, in boundary layers, we focused on a parallelizable algorithm to achieve high computational efficiency in regards of algorithmic complexity. For the presented approach, an algorithmic complexity of (N ⋅ N q ) can be estimated, where N is the number of interpolation points and N q the size of the local stencils. Since for large datasets N ≫ N q , the scaling of the presented method behaves better in comparison to classical (global) RBF interpolation, which has complexity (N 2 ). The choice of the local interpolation patches is very important for the interpolation of distorted meshes, as they occur for example in boundary layers of flow simulations. Therefore, we proposed a connectivity preserving, mesh-based neighbor search. This modified patch algorithm searches for the next directly connected neighbor elements over the common global nodes. Additionally, the neighbors of the direct neighbors are also taken into account by applying a layer strategy. To analyze the convergence of the RBF interpolation and to quantify the choice of neighbors, analytic functions were prescribed on a unit cube with different node distributions and the L 2 -error was observed after interpolating to a different mesh with equidistant node distribution. In doing so, we found that both search methods perform well for well distributed data points. However, for large aspect ratios of boundary layer meshes, the connectivity preserving, mesh-based neighbor search significantly outperformed the kd-tree based nearest neighbor search. Furthermore, constraints and boundary conditions are incorporated to enhance the interpolation of boundary layers.
Additionally to RBF interpolation, we have presented how to compute spatial derivatives based on RBFs. The derivatives were computed in a local-global approach, using a Gaussian kernel on local point stencils. Investigations on the local-local approach, using a Wendland kernel on a local point stencil, showed reduced accuracy for skewed grids. Moreover, an optimization strategy to find an optimal value of the shape parameter scaling the Gaussian kernel was presented. Good convergence of this scaled exponential kernel was shown for the computation of the gradient of an analytic field that was prescribed on a unit cube with different node distributions.
Finally, a numerical investigation with a co-rotating vortex pair was carried out to show the proposed workflow from CFD results to aeroacoustics using RBFs. This example also included the comparison to an analytic solution, which showed excellent agreement even for relatively coarse aeroacoustic meshes.