Ffowcs Williams – Hawkings analogy for near-field acoustic sources analysis

Abstract

The paper expands the scope of applying the Ffowcs Williams – Hawkings integration method. We propose using the acoustic field generated from time-dependent data stored on the FW-H control surface as the same common field for computational acoustic beamforming and dynamic mode decomposition methods to analyze the aerodynamic noise sources. We exemplify that it leads to obtaining mutually consistent and complementary information for reliable prediction of acoustic sources characteristics in the process of inverting data produced by a CFD simulation. Moreover, as the results of applying computational acoustic beamforming and dynamic mode decomposition methods depend on many geometric and algorithmic inputs, the proposed approach makes it possible to use various sets of the latter for a comprehensive analysis of obtained inversions and to form the final answer by an averaging procedure. We illustrate this by taking advantage of fast generating the examined acoustic field snapshots in any required region by the FW-H integration method for the recently developed new inverse computational acoustic beamforming algorithm and the standard dynamic mode decomposition method when carrying out a sensitivity study of the predicted acoustic source. The capabilities of the developed approach are demonstrated on the data of CFD scale-resolving simulation of turbulent flow over the 30P30N high-lift configuration.

Keywords

Beamforming computational experiment,turbulent flow aeroacoustics acoustic source dynamic mode decomposition Ffowcs Williams-Hawkings analogy 30P30N high-lift configuration

Introduction

The fundamental work of Prof. Ffowcs Willams¹ published more than a half of century ago strongly contributed to development of aeroacoustics. Further, the proposed in 1969 Ffowcs Williams – Hawkings (FW-H) analogy² provided the method for evaluation of far field acoustics in the widest range of applications in aviation industries. The FW-H integration method is based on the Lighthill’s acoustic analogy and exploits the assumption that the gas media is governed by the linear equations outside some control surface that covers the region of nonlinear acoustic source. This approach is widely used in computational aeroacoustics to analyze the distribution of sound noise in the far field.

In the paper we study the possibility of using acoustic fields calculated from time-dependent data stored on the FW-H control surface to analyze the aerodynamic noise sources by the methods of beamforming and dynamic mode decomposition (DMD). Being originally suggested for analysis of experimental data, beamforming can be also used for data sets supplied by CFD simulators, see³ and recently proposed.⁴ Such class of methods is called the computational acoustic beamforming (CAB), and we deal with CAB hereafter. Motivation for the undertaken study is as follows. First, the joint analysis of CAB and DMD results allows us for collecting mutually consistent and complementary information required for reliable prediction of acoustic sources characteristics. Second, it is well known, and we also demonstrate this, that CAB and DMD results depend on choosing geometrical parameters (placement of regions with sources, microphones, and field snapshots), discretization parameters (steps of frequency, time, and space grids), regularization and other parameters of numerical algorithms. It means that a statistical analysis of massive solution results with different CAB and DMD methods settings for a given CFD simulation is necessary to provide averaged acoustic source characteristics and uncertainty quantification. Therefore, use of the same linear acoustic field generated sufficiently fast in any required region by the FW-H integration method makes the analysis reliable and not too computationally expensive. That is why we consider this approach to be promising in computational acoustics for predicting sources of unwanted acoustic noise and here present some early results in this direction. Note that in our approach the acoustic field is considered in the near zone with respect to the FW-H control surface.

The technologies for detecting and analyzing acoustic sources are in high demand when designing modern aircraft due to constantly increasing constraints for permissible noise levels. Beamforming presents a method for localization of acoustic sources and revealing their spectral properties. To date, it is a well-established approach with a lot of theoretical papers at the background and a lot of practically useful techniques. The method of beamforming was originally developed as applied to the data of physical experiments including those performed in wind tunnels. Great efforts are made to design efficient acoustic measuring systems to properly register acoustic fields and to develop the beamforming methods to treat the corresponding data with acceptable accuracy.

The classical beamforming approaches^5,6 assume the usage of phased microphone arrays with coherent summation of the incoming signals which can be done by different methods. One of the main methods of beamforming is conventional beamforming^5,7^. It has limitations to be used in many applied problems because of the Rayleigh criterion⁸ that establishes the connection between the resolution of the source, its frequency and the geometry of the microphone array. The first ideas to address this deficiency^9,10 propose special configurations of microphone arrays. The beamforming was further developed in CLEAN¹¹ that improves the accuracy, but still has problems with the treatment of noisy signals, particularly, coherent noise. To overcome this drawback, the deconvolution techniques such as DAMAS¹² and DAMAS2¹³ were developed to refine the noise maps obtained by the conventional beamforming with the help of iterative deconvolution of the result with the so-called point spread function to take into account monopole scattering. Further improvements have resulted in more stable versions such as DAMAS-C¹⁴ and CLEAN-SC.¹⁵ Among the recent developments we would refer to functional beamforming¹⁶ that is based on using the power function on cross-spectral matrix.

Another direction of beamforming methods development is associated with the search for mathematically correct methods to solve the inverse radiation problem, in which the right-hand side of the corresponding forward boundary value problem is determined. The generalized inverse beamforming method proposed in¹⁷ is based on calculating the pseudoinverse solution of the matrix radiation transfer equation in which the matrix represents the sound radiation patterns from the source grid points to the microphone positions. The paper also considers the techniques for the regularization of the initially ill-posed system of linear equations. In the beamforming method⁴ the inverse radiation problem is formulated for the acoustic source modeled by a continuous function, as opposed to the traditional beamforming point-wise sources, and the acoustic radiation transfer operator is approximated by applying finite element techniques. Thanks to using the physically grounded parameters of grids for the source function and microphones, as well as to the high-accuracy method of calculating the integrals in the radiation transfer operator, the resulting system of discrete equations is well solvable even without regularization We use our CAB algorithm⁴ as a base for developing the approach proposed in the current paper.

The proper orthogonal decomposition (POD) and DMD methods have been actively developed since relatively recently. Their goal is the spatiotemporal analysis of a set of the flow snapshots at successive times in a certain subdomain of the flow.¹⁸ These snapshots can represent pressures, velocities and other functions resulting from a physical or numerical experiment. Both methods with several variations approximately describe a set of snapshots by a linear combination of spatial modes with time-dependent coefficients. These expansions are a convenient tool to study flow properties and, on the other hand, to reduce the dimension of the space of parameters which approximate physical fields of the examined flow with good accuracy. Particularly, the orthogonal basis used in POD¹⁹ minimizes the loss of accuracy while describing the set of considered states in it. This method reduces to the singular value decomposition which, by the Eckart – Yang theorem, gives the best low-rank approximation for the considered matrix, the snapshot matrix in our case. For example, if one uses the flow velocity as a variable, POD conserves the maximum amount of the flow energy for a given number of modes, and the modes are ranked according to the flow energy which they describe. Temporal factors for spatial modes are some general functions (mixture of frequencies) unless a special kind of method is used – the spectral POD.¹⁹

In DMD each dynamic mode is also represented by an amplitude distribution in space but without orthogonal property. The temporal factor in the spatial mode is a complex exponential with its own frequency in time and damping ratio. Hence, a linear expansion of this kind immediately gives us a certain set of frequencies (similarly to Fourier decomposition). This is one of the reasons why in the first place we consider the DMD method for analyzing the source function in the near field.

The DMD method was introduced in²⁰ for processing the results of hydrodynamic experiments and calculations. It was shown in²¹ that the dynamic modes are related to the spectrum of the Koopman operator, a linear operator describing nonlinear dynamics in the Hilbert space, which indicates that DMD can be applied to nonlinear problems as well. In²² the equivalence of DMD and the discrete Fourier transform in a certain limiting case was established as well as the method of reducing DMD dimensionality keeping the minimal possible deviation from the original data in the l₂-norm, called Optimized DMD, was described. Various modifications of the DMD algorithm were also considered in²³ to increase its robustness, in²⁴ to increase the accuracy associated with the error in data and in other works, e.g.^25,26,27

In²⁸ the currently established approach for dynamic mode calculation using the singular value decomposition was formulated, which made the method more reliable in comparison with the ill-conditioned algorithm.²⁰ It is also important that this algorithm performs a similarity transformation of the original evolution matrix, which allows to obtain its spectrum with much lower computational costs.

In the case when the original problem is strongly nonlinear, for a more accurate analysis of a dynamical system set of states, a large number of snapshots may be required for DMD, hence, the size of the evolution matrix may become too large. In this regard, the so-called kernel methods are also considered, when a certain set, called the dictionary, of transformations is selected to reduce the dimension of the basis arising in the spectrum approximation of the infinite-dimensional Koopman operator. For instance, in²⁹ a robust kernel method is presented, which allowed to separate with good accuracy the linear and nonlinear dynamics for dynamical systems of several different types.

The aim of the current paper is to develop an approach of using the CFD simulation data accumulated on FW-H surface for the subsequent analysis of sound sources in the near field by the methods of CAB and DMD.

The outline of this paper is as follows: We give a short description of the beamforming method developed in⁴ in Section Computational acoustic beamforming algorithm and present the DMD method according to³⁰ in Section Dynamic mode decomposition. Section ANALYSIS OF ACOUSTIC SOURCES GENERATED BY THE 30P30N HIGH-LIFT CONFIGURATION is devoted to the application of beamforming and DMD analysis to the data obtained in the computational experiment. As an example of the CFD simulation described in Section CFD simulation setup, we consider the turbulent flow generated around the 30P30N aerodynamic model (straight wing with high-lift configuration). Application of beamforming for the obtained data is considered in Sections Beamforming: acoustic source amplitude for the certain frequency and Beamforming: frequency-response characteristic. For the first beamforming setup, we solve the problem of finding the source function on a piecewise linear “centerline” of the wing for a certain frequency with microphones located on the FW-H surface. Then in the second setup we consider the microphones grid on an auxiliary line under the wing and look for the source function with support on a horizontal line located between the wing and microphones. The input data sets for microphones are formed by applying the FW-H integration method. To obtain the frequency-response characteristic of the source function in a certain frequency band we perform massive CAB runs. Section DMD analysis of acoustic field under the wing presents the results of the DMD analysis for the acoustic field under the same 30P30N wing. The rectangular domain with snapshots is located outside the FW-H surface, and the FW-H integration method is used to generate snapshot datasets. Several variants of DMD input data are considered in order to study the dependence of the predicted sound field characteristics on them and how close they are to those of sound sources obtained by beamforming. Finally, in the last section we discuss the results and suggest conclusions for the proposed approach of joint beamforming and DMD analysis of sound sources in near-field using FW-H analogy to generate input data.

Computational acoustic beamforming algorithm

Differential formulation of the beamforming problem

Beamforming is based on the approximation according to which the sound pressure $p$ satisfies the Cauchy problem for the wave equation in 3D space. If an object under study moves at a constant subsonic velocity, it is convenient to introduce the related coordinate system. Then the medium moves, for example, along the x₁ axis of the Cartesian coordinate system $(x_{1}, x_{2}, x_{3})$ with a Mach number $M_{\infty} < 1$ , and the equations have the form

\begin{array}{l} \frac{1}{c^{2}} \frac{\partial^{2} p}{\partial t^{2}} + 2 \frac{M_{\infty}}{c} \frac{\partial^{2} p}{\partial t \partial x_{1}} - Δ p + M_{\infty}^{2} \frac{\partial^{2} p}{\partial x_{1}^{2}} = q (x, t), x \in ℝ^{3}, t > 0, \\ p (x, 0) = \frac{\partial p}{\partial t} (x, 0) = 0, x \in ℝ^{3}, \end{array}

(1)

where a source function

q (x, t)

is nonzero on some bounded surface

S

. The surface

S

is usually located in the vicinity of the object taking into account the peculiarities of the geometry as shown in Figure 1 with example of a high-lift wing segment.

Figure 1.

The wing segment surface and the surface $S$ of acoustic source function in three-dimensional (a) and two-dimensional (b) setups.

To find the function $q (x, t)$ , an external surface $D$ with respect to $S$ is introduced, on which the values $p (x, t), x \in D$ are specified; it matches the location of the microphones.

Then the acoustic beamforming problem is formulated as the inverse problem for the right-hand side: find the source function $q (x, t)$ in equation (1) by using the values of $p |_{D}$ of acoustic field.

Assuming that harmonic oscillations are established, i.e. $p (x, t) = P (x) e^{i 2 π f t}, q (x, t) = Q (x) e^{i 2 π f t}$ , from equation (1) we obtain the equation in the frequency domain for the function P(x)

k^{2} P + Δ P - 2 i k M_{\infty} \frac{\partial P}{\partial x_{1}} - M_{\infty}^{2} \frac{\partial^{2} P}{\partial x_{1}^{2}} = - Q (x), k = \frac{2 π f}{c},

(2)

and come to the beamforming problem in the frequency domain: find the source function

Q (x)

on surface

S

in equation (2) by using the values

P |_{D}

For an arbitrary $Q (x)$ , equation (2) has the following solution³¹

\begin{array}{l} P (x) = \int_{S} Q (y) G_{M_{\infty}} (x - y) d y, x \in D, \\ G_{M_{\infty}} (x) = \frac{1}{4 π} \frac{e^{- i k^{'} (x^{'} - M_{\infty} x_{1})}}{x'}, \\ x' = \sqrt{x_{1}^{2} + (1 - M_{\infty}^{2}) (x_{2}^{2} + x_{3}^{2})}, k' = \frac{k}{1 - M_{\infty}^{2}} . \end{array}

(3)

Therefore, we use equation (3) to write the beamforming problem in the operator form

P = T Q,

(4)

by introducing the acoustic radiation transfer operator

T

, called hereafter the beamforming operator, that forms the pressure field on the surface

D

from the source function specified on the surface

S

Discrete formulation of the beamforming problem

Let ${y^{n}}_{n = 1}^{N}$ be a grid with $N$ nodes on $S$ , and $s = {(s_{1}, s_{2}, \dots, s_{N})}^{т}, s_{n} = Q (y^{n})$ a discrete vector of the source function values in the grid. Similarly, ${x^{m}}_{m = 1}^{M}$ is a grid with $M$ nodes on $D$ , which we call microphones, and $d = {(d_{1}, d_{2}, \dots, d_{M})}^{т}, d_{m} = P (x^{m})$ is the vector of the microphone signals. Usually, the inequality $M ≫ N$ is true. We denote by $T_{a} \in ℝ^{M \times N}$ a discrete approximation of the operator $T$ in equation (4); it is defined in equation (6).

The discrete beamforming problem consists of finding the vector s satisfying

{‖ \tilde{d} - T_{a} s ‖}_{2}^{2} \to min_{s},

(5)

where

\tilde{d} \approx d

is the input data (perhaps disturbed).

For approximation of operator $T$ in equation (4) we use piecewise linear interpolation of $Q (y)$ on $S$ by the values s at nodes ${y^{n}}_{n = 1}^{N}$ . Let ${ψ_{n} (y)}_{n = 1}^{N}$ be a corresponding basis elements of piecewise linear function space with values 1 at points $y^{n}$ and 0 at others. Then the source function $Q (y)$ is replaced by its piecewise linear approximation

Q (y) \approx \sum_{n = 1}^{N} s_{n} ψ_{n} (y) .

We obtain from equation (3) the formula for the components of d

d_{m} = \sum_{n = 1}^{N} s_{n} \int_{S} ψ_{n} (y) G_{M_{\infty}} (x_{m} - y) d y .

Thus, the entries of the matrix $T_{a}$ in equation (4) are calculated by the integration

T_{m n} = \int_{S} ψ_{n} (y) G_{M_{\infty}} (x_{m} - y) d y .

(6)

Note that for low frequencies such that the sound wavelength is larger than the span of the wing segment, we should exclude dependence of $Q (y)$ on the coordinate $y_{3}$ (the periodicity conditions for CFD simulation are often used for the wing segment in this direction). Therefore, it is no sense in making discretization along the coordinate $y_{3}$ for $Q (y)$ , and we consider two types of the source $S$ domain geometry: the surface and the curve, respectively, see Figure 1. The integrals in equation (6) are approximated using Gauss formulae with seven points for triangles ( $S$ is the surface) and four points for intervals ( $S$ is the curve) to achieve the required high accuracy of calculating $T_{m n}$ when the source grid is constructed in accordance with the rules of Section ‘Constructing well-conditioned beamforming matrix’.

When $M > N$ , the problem in equation (5) has a solution according to the least squares’ method

s = {(T_{a}^{*} T_{a})}^{- 1} T_{a}^{*} \tilde{d}

(7)

Constructing well-conditioned beamforming matrix

Parameters of source function and microphones grids drastically affect the condition number $C_{a}$ of the matrix $T_{a}^{*} T_{a}$ and, therefore, the error in calculating the source function according to equation (7). In this case, the fulfilment of the constraint $C_{a} < 10^{5}$ turns out to be good enough. We denote by ${dist}_{S}$ and ${dist}_{D}$ the source function and microphones grid spacing, respectively, and by ${dist}_{S D}$ some middle distance between surfaces $S$ and D. One of the main parameters of beamforming is the maximal wavelength of the examined acoustic signal. From equation (3) its value is estimated as $λ = (1 + M_{\infty}) c / f$ . Let us introduce the scaling coefficients $a_{S D}, a_{S}, a_{D}$ for the relationships

{dist}_{S D} = a_{S D} λ, {dist}_{S} = a_{S} λ, {dist}_{D} = \frac{a_{D}}{a_{S}} {dist}_{S D} .

The research results show⁴ that in the interval $a_{S D} \in [3,30]$ the following recommendations should be adhered to fulfill the constraint $C_{a} < 10^{5}$ :

• for two-dimensional case (the source function is set on a curve): $0.4 < a_{S} < 1.5, 0.08 < a_{D} < 0.15$ (e.g. we can fix $a_{S} = 0.6, a_{D} = 0.12$ );

• for three-dimensional case (the source function is set on a surface): $0.6 < a_{S} < 1.5, 0.08 < a_{D} < 0.15$ (e.g. we can fix $a_{S} = 0.8, a_{D} = 0.1$ ).

Besides, it is necessary to ensure that there are several microphones close to the node of the source grid and to satisfy the condition $M > N$ .

Note that the constraints indicated here are independent of the input data $\tilde{d}$ and are checked at the stage of generating grids and matrix $T_{a}$ .

Dynamic mode decomposition

The DMD allows to determine spatial modes in physical fields that oscillate with certain time frequencies and exponentially decay with respect to time. The problem of finding dynamic modes can be formulated as finding a linear operator of dynamic system evolution

A : p_{k + 1} \approx A p_{k}, k = 1, \dots, m,

(8)

where p is a vector of a dynamic variable in a certain set of

n

points of the analyzed domain of space; p _k is a snapshot, that is the value of p at one of the

m + 1

points in time, uniformly following each other with a time step

Δ t

. Calculation of the spectrum of the operator

A

gives the sought-for set of frequencies and the corresponding spatial modes. Let us describe the algorithm for obtaining this decomposition, according to.³⁰

First, two matrices with $n$ rows and $m$ columns are formed

X = [p_{1} | p_{2} | \dots | p_{m}] and X' = [p_{2} | p_{3} | \dots | p_{m + 1}] .

The aim of the DMD algorithm is to find the linear operator A, which approximates the evolution of considered time snapshots, that is

X' \approx A X .

(9)

Formulating this problem in the form of Frobenius norm

| | X' - A {X | |}_{F}

minimization, one can show that

A = X' X^{+}

, where X⁺ is a pseudoinverse matrix which is calculated by

X^{+} = V Σ^{- 1} U^{*}

after the singular value decomposition

A X = U Σ V^{*} .

(10)

Thus, we come to the formal expression

A = X' V Σ^{- 1} U^{*} .

Usually,

n ≫ m

. In order not to solve the spectral problem directly for the large

n \times n

matrix A, we can use the matrix

\tilde{A} = U^{*} A U,

(11)

that is the projection of A onto the orthogonal basis of columns of matrix U (onto the spatial modes of the POD method, since decomposition in equation (10) also underlies it). The operation on the right-hand side of equation (11) is called the similarity transformation, it preserves the eigenvalues.³² Therefore, the spectral problem is solved for a smaller

m \times m

matrix

\tilde{A} W = W Λ,

(12)

and the transition to the eigenvectors of A can be made by calculating the matrix

Φ = U W

, the columns of which are the required spatial modes.

Note that the computational complexity of the DMD can be further reduced by using truncated SVD.³² In this case the square matrix $\tilde{A}$ has size $r \times r$ , where $r \leq m$ is some assigned rank of the matrix in truncation operation.

Analysis of acoustic sources generated by the 30P30N high-lift configuration

CFD simulation setup

We use the 30P30N high-lift configuration to evaluate the proposed way of acoustic source analysis. The 30P30N configuration presents the three-component airfoil considered at the stage of landing. The airfoil with the slat and flap (both deployed at 30° angles) is immersed into the uniform flow at Mach number $M_{\infty} = 0.17$ with the angle of attack 5.5°. The Reynolds number based on the dimensionless airfoil chord $L = 1$ in the stowed configuration equals to $Re = 1.7 \cdot 10^{6}$ (the physical length of the airfoil chord is 0.457 m). The slat with length 0.17 has a blunt trailing edge. This setup was experimentally studied in^33,34 and taken as one of the cases at the BANC-III Workshop where the numerical results obtained by different scientific teams with the use of different codes were assessed.³⁵

For the simulations we considered a wing segment of length $0.111$ in the spanwise direction and then set the periodicity conditions. The unstructured mixed-element mesh of 34.9 M nodes was built. The predictions were carried out using the IDDES³⁶ approach implemented in the in-house code NOISEtte.^37,38 The results of simulations partly presented in³⁹ quite well correspond to the available experimental data. Figure 2(a,b) show the instantaneous field of pressure time derivative (in black and white shades) that reveals the presence of acoustic sources with sound waves of different lengths which are localized at different parts of the configuration: in particular, around the slat-wing and wing-flap gaps, and the trailing edge of the flap. Besides, one can see the turbulent content colored in velocity magnitude values and the FW-H surface where all the needed data is accumulated. Figure 2(c,d) show the computed and experimental pressure spectrum measured at the points P2 and P6 located on the slat surface (we use the non-dimensionalized notation: the Strouhal number $S t = f L / U_{\infty}$ , the frequency f, the upstream velocity $U_{\infty} = 1$ , the dynamic pressure $q_{\infty} = ρ_{\infty} U_{\infty}^{2}$ , the ambient density $ρ_{\infty} = 1$ , the sound speed $c = 1 / M_{\infty}$ ). In both pictures one can notice peaks roughly corresponding to the Strouhal numbers St = 11, St = 15, and St = 21. These peaks are clearly seen in the experimental data and quite well, up to a slight shift to the right, caught by the numerical prediction.

Figure 2.

30P30N case simulation results. Instantaneous field of the pressure derivative with respect to time are in gray, vorticity structures are colored by the velocity magnitude, the FW-H surface is drawn in yellow (a); A zoom of the slat region with the indicated points P2 and P6 of hotwire measurements (b); Power spectral density computed in the points P2 and P6 (c and d, respectively).

Beamforming: acoustic source amplitude for the certain frequency

Consider first the analysis results for one frequency, $S t = 21$ . We use a piecewise flat surface $S$ to define the acoustic source function on it, Figure 3(a). The microphone grid is placed on the FW-H surface. The computational experiment data accumulated on the FW-H surface is averaged over time windows and filtered in the one-third octave frequency band 18.7 < St < 23. As the wavelength $λ = M_{\infty}^{- 1} S t_{c}^{- 1} \approx 0.28$ for the central frequency St_c = 21 exceeds the span length 0.111 of the wing segment, the source function must be constant along the x₃ axis. This means that beamforming can (and should) be done using 2D setup in this example, cf. Figure 1. We actually tested both 2D and 3D setups, similarly as it was done in.⁴ The 2D and 3D beamforming results shown in Figure 3(b) are consistent with each other and predict the presence of a sufficiently strong acoustic source in the region $x_{1}$ ≈ 0.1 which corresponds to the slat-wing gap.

Figure 3.

The microphone surface D and the piecewise flat source surface S with the source amplitude distribution for 3D setup (a); amplitude of the source function on the polyline for the 2D setup and along the centerline of the piecewise flat surface for the 3D setup (b).

Beamforming: frequency-response characteristic

Now let us consider the application of the developed beamforming algorithm to evaluate the frequency-response characteristic (FRC) of the acoustic source in the desired frequency range, in this example, the range 8 < St < 35. The lower part of half-space below the wing (sound radiation towards the Earth’s surface) is of interest. Therefore, for 2D beamforming setup, the microphones curve D (half-ellipse) and source function line S ( $- 0.2 \leq x_{1} \leq 1.7$ , $x_{2} = - 0.1$ ) are selected in the way shown in Figure 4. The data for the microphone array at $D$ is generated by the FW-H integration method, followed by averaging over time windows and filtering in the considered frequency band.

Figure 4.

Acoustic source function line S and microphones array curve $D$ for FRC analysis.

To obtain the FRC, we perform massive beamforming on the sampled frequencies generated by subdividing the band 8 < St < 35 evenly on a logarithmic scale in 1/12-octave increments. Note that for each sampled frequency the own source function grid is used in accordance with the rules of Section Constructing well-conditioned beamforming matrix. Moreover, beamforming for each frequency consists of several calculations using grids for close frequencies and then averaging the piecewise linear functions thus obtained. Figure 5(a) shows these averaged plots of the source functions for three sampled frequencies St = 11.4, 15.3, 21.5 as example. The rise in the amplitude of the acoustic source for the coordinate corresponding to the region of slat-wing gap is clearly noticeable in all plots. Also, some changes in the amplitude are observed behind the wing. The FRC in Figure 5(b) is estimated by calculating the maximum amplitude of the source function at sampled frequencies. Note that the observed maxima of FRC about St = 11, 15.5, 21.5, 27 are caused by maxima of the source function plots for these frequencies in the region of the slat-wing gap. The resulted FRC is consistent with the pressure spectrum predicted at a point near the slat-wing gap, see Figure 2(c,d), showing peaks at the same frequencies.

Figure 5.

Acoustic source function. Plots for frequencies St = 11.4, 15.3, 21.5 (a); frequency-response characteristic with the noise maxima labeled by vertical lines at St = 11,15.5, 21.5, 27 (b).

DMD analysis of acoustic field under the wing

The solution of the beamforming problem allowed us to obtain the FRC of the acoustic pressure from the wing and approximately estimate the location of the sources along the chord. Distinct noise maxima observed in the regions St = 11, 15.5, 21.5, 27 are obviously of primary interest in studying the acoustic noise generated by the flow. Here we explore the possibility of analyzing the properties of the sound source in the slat-wing gap using DMD method and show that these results can successfully complement the information obtained from the beamforming.

First, note that, although the problem is three-dimensional, it is sufficient to consider a certain 2D domain with snapshots in the plane of the wing cross-section, i.e. not a 3D one. This simplification of DMD analysis for large acoustic wavelengths compared to the wing segment length is similar to the application of beamforming in the 2D setup, see Section ‘Discrete formulation of the beamforming problem’.

Each snapshot $p_{k}$ is a vector of values of the sound pressure p, evaluated by using the FW-H integration method at grid points of the flow region under study. We consider two variants for the region with snapshots: rectangles P1 and P2 located closer to the wing and farther from it, Figure 6.

Figure 6.

Position of the airfoil, FW-H surface and rectangular region P₁ $(x_{1} \in [- 0.3, 1.2], x_{2} \in [- 0.75, - 0.25])$ with a grid on which snapshots are calculated. Red border also shows the shifted down region P₂ $(x_{1} \in [- 0.3, 1.2], x_{2} \in [- 1.0, - 0.5])$ .

The grids in P1 and P2 have $n \sim 10^{4}$ points for a good resolution. The number of snapshots $m = T / Δ t$ for the time interval $T$ is estimated as follows. We take $T$ equal to about 10 periods of the minimum frequency from the range of interest $8 < St < 35$ and take about 50 snapshots for the maximum frequency. Thus, we come to the estimate $m \sim (10 / S t_{min}) \times (50 S t_{max}) \sim 2000$ .

After forming matrices $X, X'$ , obtaining the spatial modes $Φ$ and the corresponding eigenvalues $Λ$ for the operator $A$ from equation (9), we can analyze the dynamics of the spatial modes which characterize the sound pressure distribution on the grid in the rectangle P₁ or P₂.

According to equations (8), (11), and (12) the vector of sound pressure in the considered rectangle at an arbitrary time $t$ can be approximated by the expression

p (t, x) \approx \sum_{j = 1}^{m} P_{j} φ_{j} (x) e^{(σ_{j} + i ω_{j}) t},

(13)

where

φ_{j}

is a spatial mode from

Φ

with some coefficient

P_{j} \in ℂ

in the expansion. The frequencies

ω_{j}

and damping ratios

σ_{j}

are found by the formula

σ_{j} + i ω_{j} = \log (Λ_{j}) / Δ t

We use a 150 × 80 grid in each of the rectangles P₁, P₂ for snapshots (i.e. n = 12000) and $Δ t_{0} = 8 \cdot 10^{- 4}$ as the minimal time interval between them. First, we will discuss the issue of the reliability of the DMD analysis results, i.e. the dependence of the dynamic modes on the method for generating snapshots.

Figure 7 shows the plots of the singular values in descending order for the decomposition in equation (10) in two regions P₁ and P₂. One can see their distinct similarity (here $Δ t = Δ t_{0}, m = 2400$ ).

Figure 7.

Singular values of the pressure field snapshot matrix $X$ .

The DMD method, unlike the Proper Orthogonal Decomposition, does not allow to determine which of the modes calculated to describe the acoustic field are more significant than others. However, we can select those modes that damp most slowly. Figure 8 displays frequency values $f_{j} \equiv St = ω_{j} / 2 π$ and logarithmic decrements $δ_{j} = σ_{j} / f_{j}$ in the interval $16 < f_{j} < 29$ and $0 > δ_{j} > - 0.02$ respectively, i.e. values for the weakly damped modes obtained in the regions P₁ and P₂.

Figure 8.

Distribution of points $(f_{j}, δ_{j})$ for weakly damped modes in P₁ and P₂.

We see that for both regions there are frequencies near $St = 17, 21, 27$ (vertical lines) which correspond to the set of high-amplitude frequencies $St = 15.5, 21.5, 27$ obtained by beamforming, Figure 5(b). Thus, complementing the solution of the beamforming problem with DMD results allows to more accurately identify physically important modes and then analyze their distribution in space.

However, it is also important to check that dynamic modes with weak damping make a sufficient contribution to the total sum in equation (13) for their interpretation to have sense. Figure 9 shows the relative weights of these modes (marked with triangles) in the dynamic mode decomposition for the frequency interval $14 < St < 29$ ; the weights represent the ratio of modulus of amplitude to the sum of modules of all amplitude from the interval. Their contribution is quite noticeable considering that the other modes marked with circles are characterized by strong damping.

Figure 9.

Relative amplitude weight for all modes in the interval 14 < St < 29. Triangles represent the weakly damped modes from Figure 8.

To estimate the sensitivity of DMD analysis to the parameter settings, we changed the sets of snapshots. Figure 10 shows DMD analysis results in P₁ for two other time steps while generating snapshots: with $Δ t = 2 Δ t_{0}$ and $Δ t = 3 Δ t_{0}$ .

Figure 10.

Weakly damped modes for snapshot sets in P₁ with parameters $Δ t = 2 Δ t_{0}, m = 2400$ (blue triangles) and $Δ t = 3 Δ t_{0}, m = 1600$ (red triangles).

We see from Figures 8, and 10 that, although the frequencies and decrements obtained differ slightly in each of the four parameter settings, there are distinct weakly damped modes with frequencies close to $St = 17, 21, 27$ . Moreover, the results with $Δ t = 2 Δ t_{0}$ and $Δ t = 3 Δ t_{0}$ are shifted to the left relative to the noted frequency $St = 17$ and they are closer to the beamforming peak frequency St = 15.5, Figure 5(b). One way or another, the observed variability suggests that it is necessary to carry out several DMD calculations on different sets of snapshots and perform the appropriate statistical analysis.

Figures 11–13 show six spatial modes with frequencies near St = 17, 21, 27 corresponding to the weakly damped modes $(f_{j}, δ_{j})$ in Figure 8 for P₁ and P₂.

Figure 11.

Dynamic mode with the frequency St = 17.4 in Р₁ (a) and St = 17.2 in Р₂ (b).

Figure 12.

Dynamic mode with the frequency St = 21.1 in Р₁ (a) and St = 20.7 in Р₂ (b).

Figure 13.

Dynamic mode with the frequency St = 26.9 in Р₁ (a) and St = 26.8 in Р₂ (b).

We may conclude that it is preferable to locate the region for DMD acoustic source analyses as close to the airfoil as possible, i.e. use Р₁ in our case, to reveal the fine structure of the source function.

Note also that the amplitude distribution in the Р₁ region suggests that we are mainly dealing with a monopole-type source located at the bottom near the slat-wing gap. To confirm this consideration, we performed an additional beamforming calculation at the frequency St = 17.4 and estimated the acoustic field in Р₁ for the found source function. The similarity of both patterns can be seen in Figure 14.

Figure 14.

Acoustic field distribution in P₁ at the frequency St = 17.4 corresponding to the DMD mode (a) and to the source function predicted by CAB (b).

Discussion and conclusions

We propose the approach to analyze distributed sound sources in the flow near field using jointly the method of CAB and DMD, both methods exploit the FW-H analogy to provide the same input data. This approach is based at least on two following considerations. Firstly, the FW-H analogy itself assumes that the linear acoustic model is valid outside the FW-H surface. So, both the CAB and DMD methods that are aimed at obtaining exact solutions of linear problems (wave equation and linear dynamical system correspondingly) can provide consistent results. Secondly, both methods provide the information about the sound field and acoustic sources in mutually complementary forms, which allows to improve the prediction accuracy. Note that the CAB and DMD algorithms are sensitive to the noise of the FW-H surface data caused both by the “nonlinearity intrusion” and numerical errors. This unwanted noise is mitigated by averaging results of massive runs of these algorithms with different settings.

To evaluate the applicability of the developed approach, we use the numerical results of a 3D CFD simulation of the unsteady turbulent flow over a 30P30N high-lift three-component configuration of the aircraft wing during landing. The numerical data was previously validated against the experimental results and reference numerical results of other authors. The corresponding unsteady data was stored at the FW-H surface and then was transferred by the FW-H integration method onto grids of virtual microphone arrays within the CAB procedure,⁴ and onto the snapshots grids when applying the DMD according to the algorithm.³⁰ Since the span length of the 30P30N segment considered in the CFD simulation was less than the sound wavelength from the frequency band of interest $8 < St < 35$ , we used 2D setups for the analysis, namely, source function curves in beamforming, and plane snapshots in DMD.

When doing the beamforming analysis, we consider two geometric setup variants with different source function curves: “inside” and under the airfoil, respectively. The curves for microphone arrays are also different: the FW-H surface and a half ellipse, respectively. There is a good agreement between the plots for both setups, Figures 3 and5(a), of the predicted acoustic source function with a distinct rise in amplitude near the coordinate corresponding to the slat-wing (the plots do not have to be same, since the source curves are different), which allows us to conclude the quite acceptable accuracy of the results obtained. The goal of the second setup is to calculate the frequency response characteristic of the frequency band $8 < St < 35$ (unlike the first one considered for the single frequency St = 21). To do this, we perform massive CAB runs by dividing the band into approximately 25 sampled frequencies. The predicted FRC has distinct maxima near St = 11, 15.5, 21.5, 27, see Figure 5(b); all of them are caused by the maxima of the acoustic source function plots in the region of slat-wing gap.

While performing DMD, we first looked at how stable the analysis results are when varying the data coming to the input of the algorithm. The four considered parameter settings for generating snapshot sets result in close frequency characteristics of the calculated sound fields (only low attenuation modes are compared). The best match with the maxima of the FRC of the beamforming is provided by performing DMD in the P₁ region, and with the rather small time step to form a set of snapshots. Therefore, we come to the conclusion that it is preferable to place the DMD region for source analysis as close to the wing as possible (to the FW-H surface). Analysis of the amplitude maps of spatial modes allows us to conclude that the monopole plays the main role in the source function, Figures 11–14, which justifies the use of the beamforming method setup with the monopole source function.

The results of the undertaken research show that the proposed approach to joint analysis of sound sources by applying near-field CAB and DMD methods, based on using time-dependent input data stored on the FW-H surface, can accurately predict the space-frequency characteristics of an acoustic noise source. The stability of the prediction and the corresponding uncertainty quantification require statistical analysis of the massive solution results with different settings in CAB and DMD algorithms.

Footnotes

Acknowledgements

We thank our colleagues Alexey Duben and Pavel Rodionov for providing the data of the CFD simulation and the pictures used in the paper.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Russian Foundation for Basic Research (Project 19-51-80001). The computations were carried out using the supercomputers K60 and K10 of Collective Usage Centre of the Keldysh Institute of Applied Mathematics of the Russian Academy of Sciences.

ORCID iD

Gleb Plaksin

References

Ffowcs Williams

. Noise from convected turbulence. J Acoust Soc Am 1961; 33(11): 1675–1675.

Ffowcs Williams

Hawkings

. Sound generation by turbulence and surfaces in arbitrary motion. Philos Trans A Math Phys Sci 1969; 264(1151): 321–342.

Lien

Yee

. Computational acoustic beamforming for noise source identification for small wind turbines. Int Sch Res Notices 2017; 2017: 7061391.

Kozubskaya

Plaksin

Sofronov

. Statement of the beamforming problem and a method of its solution for the localization of an acoustic source based on computational experiment data. Comput Math Math Phys 2021; 61(11): 1864–1885.

Johnson

Dudgeon

. Array signal processing: concepts and techniques. London: Prentice-Hall, 1993, p. 533.

Dougherty

. Beamforming in acoustic testing. In: Mueller

(ed) Aeroacoustic measurements. Berlin: Springer-Verlag, 2002, pp. 62–97.

Sijtsma

. Acoustic beamforming for the ranking of aircraft noise. In: Schram

Dénos

Lecomte

(eds) Accurate and efficient aeroacoustic prediction approaches for airframe noise, VKI Lecture Series 2013-03, Von Karman Institute for Fluid Dynamics. Rhode-Saint-Genèse, BE, 2013.

Carroll

Ostlie

. An introduction to modern astrophysics. International ed of 2nd revised ed. Harlow, UK: Pearson Education Limited, 2013.

Sarradj

. A generic approach to synthesize optimal array microphone arrangements. In: 6th Berlin beamforming conference. Berlin, Germany, 29 February–1 March 2016. paper no. BeBeC-2016-S4.

10.

Malgoezar

AMN

Snellen

Sijtsma

, et al. Improving beamforming by optimization of acoustic array microphone positions. In: 6th Berlin beamforming conference. Berlin, Germany, 29 February – 1 March 2016. paper no. BeBeC-2016-S5.

11.

Högbom

. Aperture synthesis with a non-regular distribution of interferometer baselines. Astron Astrophys Suppl Ser 1974; 15: 417–426.

12.

Brooks

Humphreys

. A deconvolution approach for the mapping of acoustic sources (DAMAS) determined from phased microphone arrays. J Sound Vib 2006; 294: 856–879.

13.

Dougherty

. Extension of DAMAS and benefits and limitations of deconvolution in beamforming. In: 11th AIAA/CEAS Aeroacoustics Conference. Monterey, CA, USA, 23–25 May 2005. paper no. AIAA-2005-2961.

14.

Brooks

Humphreys

. Extension of DAMAS phased array processing for spatial coherence determination(DAMAS-C). In: 12th AIAA/CEAS Aeroacoustics Conference (27th AIAA Aeroacoustics Conference). Cambridge, MA, USA, 8–10 May 2006. paper no. AIAA-2006-2654.

15.

Sijtsma

. CLEAN based on spatial source coherence. Int J Aeroacoust 2007; 6(4): 357–374.

16.

Dougherty

. Functional beamforming for aeroacoustic source distributions. In: 20th AIAA/CEAS aeroacoustics conference. Atlanta, GA, USA, 16–20 June 2014. paper no. AIAA-2014-3066.

17.

Suzuki

. Generalized inverse beam-forming algorithm resolving coherent/incoherent, distributed and multipole sources. J Sound Vib 2011; 330: 5835–5851.

18.

Taira

Brunton

Dawson

, et al. Modal analysis of fluid flows: an overview. AIAA J 2017; 55(12): 4013–4041.

19.

Holmes

Lumley

Berkooz

, et al. Turbulence, coherent structures, dynamical systems and symmetry. 2nd ed. Cambridge University Press, 2011.

20.

Schmid

. Dynamic mode decomposition of numerical and experimental data. J Fluid Mech 2010; 656: 5–28.

21.

Rowley

Mezic

Bagheri

, et al. Spectral analysis of nonlinear flows. J Fluid Mech 2009; 641: 115–127.

22.

Chen

Rowley

. Variants of dynamic mode decomposition: boundary condition, Koopman, and Fourier analyses. J Nonlinear Sci 2012; 22: 887–915.

23.

Scherl

Strom

Shang

, et al. Robust principal component analysis for particle image velocimetry, 2019.

24.

Dawson

Hemati

Williams

, et al. Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition. Exp Fluids 2016; 57: 42.

25.

Wynn

Pearson

Ganapathisubramani

, et al. Optimal mode decomposition for unsteady flows. J Fluid Mech 2013; 733: 473–503.

26.

Jovanovic

Schmid

Nichols

. Sparsity-promoting dynamic mode decomposition. J Phys Fluids 2014; 26(2): 024103.

27.

Erichson

Brunton

Kutz

. Compressed dynamic mode decomposition for real-time object detection. J Real-Time Image Proc 2016; 16: 1479–1492.

28.

Schmid

. Dynamic mode decomposition of numerical and experimental data. J Fluid Mech 2010; 656: 5–28.

29.

Baddoo

Herrmann

McKeon

, et al. Kernel learning for Robust dynamic mode decomposition: linear and nonlinear disambiguation optimization. (LANDO), 2021.

30.

Rowley

Luchtenburg

, et al. On dynamic mode decomposition: theory and applications. J Comput Dyn 2014; 1(2): 391–421.

31.

Lee

. A direct boundary integral formulation for acoustic radiation in a subsonic uniform flow. J Sound Vib 1994; 175: 51–63.

32.

Trefethen

Bau

III . Numerical linear algebra. SIAM, 1997, p. 373.

33.

Pascioni

Cattafesta

Choudhari

. An experimental investigation of the 30P30N multi-element high-lift airfoil. In: 20th AIAA/CEAS Aeroacoustics Conference. Atlanta, GA, USA, 16–20 June 2014. paper no. AIAA-2014-3062.

34.

Pascioni

Cattafesta

. Unsteady characteristics of a slat-cove flow field. Phys Rev Fluids 2018; 3(3): 034607.

35.

Choudhari

Lockard

. Assessment of slat noise predictions for 30P30N high-lift configuration from BANC-III workshop. In: 21st AIAA/CEAS Aeroacoustics Conference. Dallas, TX, USA, 22–26 June 2015. paper no. AIAA-2015-2844.

36.

Spalart

Deck

Shur

, et al. A new version of detached-Eddy simulation, resistant to ambiguous grid densities. Theor Comp Fluid Dyn 2006; 20(3): 181–195.

37.

Abalakin

Bakhvalov

Gorobets

, et al. Parallel research code NOISEtte for large-scale computations of aerodynamics and aeroacoustics problems. Vychisl Metody Programm 2012; 13(3): 110–125. (in Russian).

38.

Gorobets

Bakhvalov

. Heterogeneous CPU+GPU parallelization for high-accuracy scale-resolving simulations of compressible turbulent flows on hybrid supercomputers. Comput Phys Commun 2022; 271: 108231.

39.

Abalakin

Duben

Gorobets

, et al. Numerical simulation of slat noise using immersed boundary method on unstructured meshes. In: 25th AIAA/CEAS Aeroacoustics Conference. Delft, The Netherlands, 20–23 May 2019. paper no. AIAA-2019-2461.