Sound shielding simulation by coupled discontinuous Galerkin and fast boundary element methods

ABSTRACT A code coupling has been established for performing efficient fan tone shielding simulations of aerial vehicles with unconventional engine installations. In particular, the Fast Multipole Boundary Element Method (FM-BEM) which is formulated to solve a surface integral based on the Kirchhoff-Helmholtz wave equation for large geometries is combined with a volume resolving Discontinuous Galerkin (DG) method which is well suited for the compact region around a jet engine intake where strong mean flow gradients are present. The Möhring-Howe acoustic analogy is utilised during the backward data exchange process for derivation of acoustic velocities in presence of a mean flow. The method can help to overcome a major difficulty related to computational complexity when solving fan noise shilding and scattering problems for a complete aircraft geometry.


Introduction
Noise reduction comes into play alongside electromagnetic scattering and signature detection when dealing with problems of unmanned aerial vehicles (UAV). A multi-level fast multipole boundary element method (subsequently referred to as FM-BEM) has been implemented into the FMCAS code of DLR. The code can be used both in thread parallel and MPI modes and has an outstanding performance when it comes to solving large scale aircraft scattering and shielding problems (Lummer et al., 2013). For example, a fast BEM applied by Papamoschou (2010) and Papamoschou and Mayoral (2013) has proved to be valuable for computing jet source acoustic diffraction. Recently, Hu et al. (2017), Nark and Hu (2021) made considerable progress in developing a time domain fast BEM method. Also, for problems requiring solutions in the frequency domain, fast BEM can be very effective on its own, as demonstrated in the work of Wolf and Lele (2011), providing that sound sources are well defined and mean flow effects are of minor importance. However, because of the aforementioned limitations, the FM-BEM cannot be used for solving installed fan noise problems. In that case, a volume resolving Discontinuous Galerkin (DG) method may be well suited for computing the wave propagation through a highly non-uniform flow with a wide range of Mach numbers, typically found in the vicinity of a turbofan engine intake. With that being said, applying CONTACT S. Proskurov stanislav.proskurov@dlr.de volume resolving methods to a full-scale aircraft requires substantial high-performance resources as of today. A high level of detail and the requirement to solve for frequencies up to several kHz result in a soaring effort, making such approach unsuitable for design. The computational effort can be reduced to some extent by using an integration surface and applying the Ffowcs Williams and Hawkings (1969) (FWH) integral method for determining the far-field pressure. For example, see the work of Ren et al. (2022) where FEM was combined with acoustic analogy for simulating a complicated centrifugal fan problem. In this context, a coupled DG / FM-BEM can be regarded as an extension of a deterministic integral method where it is easily possible to enclose all diffracting airframe surfaces.
In the field of Computational Aero-Acoustic (CAA) the use of coupled methods is not a novel concept. For example, Casenave et al. (2014) applied the BEM to FEM coupling to solve an acoustic cylinder flow problem. Liu et al. (2021) investigated acoustic-structure interactions in combustion chambers featuring rigid/elastic baffles with the dual-reciprocity BEM coupled to FEM, where the latter was used to model the structural dynamics of elastic walls. In our work, the DG method is coupled to FM-BEM where a fast multipole formulation is better suited for solving engineering problems. The coupled method should be able to efficiently overcome the problem of high computational complexity when there is a significant disparity between intake and free-stream flow. In addition to that, for unconventional engine installations it can be challenging to isolate the near-field. For that reason an efficient implementation of a two-way strong coupling is presented which makes our approach unique. The simulation begins in a DG sub-domain where the Acoustic Perturbation Equations (APE) (Ewert & Schröder, 2003) are solved on a RANS background flow using a relatively small bounded domain. Thus, the global boundary conditions can be enforced on the DG subdomain via feedback from the FM-BEM. Here, the fluctuating enthalpy is used as a mediator between both worlds which is wrapped up in the Möhring-Howe acoustic analogy (Heitmann et al., 2016). During the feedback phase, the Fast Multipole Method (FMM) is utilised for efficient data plotting on all degrees of freedom at permeable DG boundaries. Our investigation of different engine placements indicate that feedback can be important if strong reflections are produced by some neighbouring aircraft components. At the end of a full cycle, the updated nearfield will have the new source information available for the FM-BEM. Then, a more accurate global solution can be established via an iterative process where the two methods complement each other by exploiting the full potential in terms of performance and fidelity. This was demonstrated at the 2021 Aviation Forum by applying the coupled method for solving a full-scale flow-acoustics problem (Proskurov et al., 2021).

Boundary integral equation
The FM-BEM simulations are performed using DLR's FMCAS code (Fast Multipole Code for Acoustic Shielding) (Lummer, 2019). This section introduces the reader to the fundamental physical concept of the method. In technical acoustics, sound waves can be regarded as small disturbances generated by a turbulent flow and if viscous effects are neglected, the relationship between pressure and velocity can be described by the Euler equation. In CAA, usually, the linearised form (LEE) is applied. However, the LEE may also contain nonacoustic disturbances, and thus, special vorticity-free version called the Acoustic Perturbation Equations (APE) can be used instead. The APE equations were derived by Ewert and Schröder (2003) to eliminate the problem of unbounded growth of hydrodynamic instabilities and in contrast to the LEE rely on the isentropic assumption. This is generally valid for computing the acoustic wave propagation both at the near-and far-field. For the FM-BEM, the APE equations have been reduced to the classical Helmholtz equation via the Taylor transformation (Taylor, 1978). Assuming the Mach number is small compared to 1, a divergence-free mean flow implies ∇ · v 0 = ∇ 2 0 = 0 where v 0 = ∇ 0 and v = ∇φ are the definitions of a mean and acoustic potential. Using the above quantities on the L.H.S. of Equation (1), namely the convective part of the APE-4 equations, as well as applying a simplified 1 convection operator D t ≡ ∂ ∂t + ∇ 0 · ∇ and substituting for pressure, yields the convective wave equation for the velocity potential: In the next step, Equation (2) is divided throughout by c 2 0 and expanded into the governing acoustics equation for steady isentropic potential flow at low Mach number, where c 0 is the speed of sound at infinity, t the time, M ≡ v ∞ c 0 is the Mach number at infinity, ≡ 0 v ∞ is the mean potential divided by the speed of the uniform stream at infinity, ( 0 ≡ c 0 M ) and φ is the acoustic velocity potential which is subject to the added source S . In the derivation of Equation (3), the term (∇ 0 · ∇) 2 φ has been omitted due to an assumption of a slowly varying mean gradient (on acoustic scale) which is applied together with the Mach number restriction. Following the transformation of Taylor (1978), the ansatz function for the acoustic potential reads: φ (x, t) = ψ (x, ω)e iωt e iκM (x) whereψ is the Taylor transformed acoustic velocity potential, κ = ω c 0 is the wavenumber and i − the imaginary unit. The source term becomes S (x, t) =Ŝ (x, ω)e iωt . The acoustic sources must satisfy a causality condition which ensures the behaviour of an outgoing wave at infinity. Substituting these into Equation (3) results in the Taylor transformed convected wave equation, where T (x, ω) ≡ e iκM 1 (x) stands for the Taylor phase factor and x ≡ (x, y, z) are the space coordinates. The following transformation rules apply: In the second equation of Equation (5), the last term is dropped by assumption, |∇ | 2 = 0. Then, by using the properties of Equations (5) in (4) and dropping O(M) 2 terms once again, we arrive at the Helmholtz equation for the Taylor transformed acoustic velocity potential ψ (x, ω), where + is the exterior CAA domain. Consider a free field Green's function, G(x, y) = e −iκr 4π r where r 2 = (x − y) 2 denotes the separation distance between two points in space (Agarwal & Morris, 2006). The Helmholtz equation for the Green's function is (κ 2 + ∇ 2 x )G(x, y) = −δ(x − y). The boundary integral equation can be obtained following the standard procedure (Crighton et al., 1992).
In Equation (7), ξ (x) is a free term coefficient which depends on the location of point x ({0} outside the domain, {1/2} on a surface or {1} anywhere in a volume inside the domain) andψ V is the incident potential term which accounts for the source. In a classical BEM approach, the surface velocity potentialψ (y, ω) is derived from an incident source which is present somewhere in the volume. In addition to that, it is possible to directly prescribeψ and ∂ψ ∂n on some part of a scattering body, denoted as d s (e.g. via interpolation from DG). Then, d s can be treated as a source surface which can entirely substitute the necessity of having theψ V term. Hence,ψ (x, ω) can be obtained from a partially known surface solution. Also, at this stage it is convenient to relate the gradient on the surface ∂ψ (y,ω) ∂n y toψ (y, ω) by introducing the proportionality function Y(y), The specific wall admittance can be related to a complex reflection factor which is the amplitude ratio between reflected and incident waves. The proportionality function, Y, is zero at solid walls whereas anywhere in free space it will be a complex number. Using the above property we may apply the differential operator D y ≡ ∂ ∂n y − Y(y) to Equation (7) which then yields for a point x on the boundary: Notice that in Equation (9) the R.H.S. is comprised of two types of sources, a volume source followed by a surface source distributed over a surface d s , with the latter not following directly from Equation (7) but derived analogously to the L.H.S. surface integral. In theory, Equation (9) is the governing BEM equation which can be discretised and solved. However, for most exterior boundary value problems the solution of Equation (9) is not unique if the wavenumber κ is an eigenvalue of the inner problem of the body with surface ∂ . One possible solution to overcome the uniqueness problem is to consider the Burton-Miller approach (Burton & Miller, 1971) where a linear combination of Equation (9) and its normal derivative w.r.t. x is used. This is equivalent to applying a simple operator to the above equation, B x = 1 + α ∂ ∂n x where a common choice for a non-zero coupling constant, α ≈ i/κ.

Discretisation
The boundary integral equation, Equation (10), has to be reduced to a system of linear equations for the surface valuesψ (x, ω) in order to be solvable numerically. First, a surface is discretised into N triangles: Second, while solving for one ω at a time, a constant value ψ j can be assumed on each element. This allows taking it outside of the boundary integral in Equation (10) simplifying the summation. Then, using the quadrature formula with collocation points y jm and integration weights W m j on ∂ j , the summation operator is defined as follows: The system of linear equations for ψ j reads: where C x ≡ 1/2[1 + αY(x)] uses the properties of wall admittance combined with the Burton-Miller operator.
The above system is solved with the FMM procedure using the plane wave expansion for the free field Green's function in Equation (12). The far-field signature is related to the near-field signature via a function which is known as the transfer function (Lummer, 2019).

Implementation of the FMM
The FM-BEM represents a fast summation method which is applied to solve a linear system in Equation (13) reducing the number of matrix-vector product operations from O(N 2 ) which is equivalent to solving an 'N-body problem' to O(N 3/2 ) for a single-level FMM and eventually to O(N log N) for a multi-level FMM, where N denotes the number of surface elements. For the FMM implementation, a domain has to be subdivided into a set of cubes. If the solution is to be determined for any chosen cube the basic idea is to perform a matrix-vector product with all its neighbours and use a suitable far-field approximation for the remaining cubes. In doing so, the computational complexity can be reduced from O(N 2 ) to O(N 3/2 ) where the FMM uses transfer of the signature functions on cubes of equal size. For the multi-level FMM (MLFMM) the scattering geometry is placed inside large cubes which are recursively sub-divided into child cubes in the form of an octree. The far-field signatures are moved up the octree from the child cubes to the centre of the parent cube. The accumulated information has to be carefully interpolated onto the collocation points of the parent cube (Orszag, 1974). The near-field signatures are transported in the opposite direction, being passed down to the child cubes from the top level. This results in a very efficient algorithm ≈ O(N log N). (See Figure 1) Many details and methodologies have to be taken into account for an efficient numerical implementation of the MLFMM which cannot be covered in this paper. Since the FMCAS code is almost a standard implementation it is possible to refer to the literature (Sylvand, 2003;Wolf & Lele, 2009, 2010. In our version of the MLFMM the iterative solvers from the Portable Extensible Toolkit for Scientific Computation (PETSc) library are used for improved performance and code optimisation. In particular, the GMRES solver was preferred for its superior convergence property that was employed together with the Generalised Additive Schwarz Method (GASM) (Goossens et al., 1996) to precondition the system of linear equations on a triangulated grid.

Discontinuous Galerkin method
The following section outlines general properties of the DISCO++ DG code (Mößner et al., 2018), mainly focussing on the features relevant to DG / FM-BEM coupling. Fundamentally, DG is a high-order method used for solving partial differential equations on unstructured meshes. The DG method used in this work, which solves the APE system (Ewert & Schröder, 2003) in a steady nonuniform mean flow, has been converted into the frequency domain to cope with the frequency domain formulation of the FM-BEM. The additional term in the governing Equation (14) that includes iω factor is a result of time to frequency space conversion where a conven- In Equation (14),q = {pûvŵ} T is a complex variable and A, B and C are flux matrices. Multiplication of the APE system with a test function, integration over volume and application of the Gauss theorem yields a weak formulation of the APE equations as in the standard DG approach. A quadrature-free implementation is achieved via linear mapping of volumes onto tetrahedral cells. There are ten degrees of freedom per face (four nodes on each edge and one face-central) which make up a total of 20 collocation points per tetrahedral cell. This allows for using third-degree polynomials with nodal shape functions (Mößner et al., 2018). In the current implementation, the method is used without padaptivity where a desired frequency resolution governs the refinement and coarsening of the mesh. Since the frequency solver only deals with a single frequency at a time, such compromise is well justified, also, improving the usability of the DG method while expanding its range to new research problems, solution advancing schemes and equations. The coupling concept presented in this paper, once implemented, applies to DG methods of arbitrary order of accuracy in space.
After discretising the governing equation with the DG method, the linear system of equations is solved with a pseudo-time marching method. Integrating with the 4th order Runge-Kutta scheme (RK4) robustly solved all cases tried up to date. As a matter of fact, the time marching accuracy is irrelevant for frequency solutions where stability and convergence are the decisive factors. Therefore, the solution path followed by means of either explicit or implicit schemes will not affect the quality of the final result providing the same convergence criteria have been reached.

Boundary conditions
The DISCO++ code features several flux reconstruction schemes which can be used for determining the DG solution between unstructured elements. The default choice for the solution of APE equations is to compute the upwind flux via characteristic splitting. First, the left and right Eigenvector matrices are defined in the rotated local coordinate system such that the normal vectorn points from left to right. Second, depending on the normal velocity the correct set of matrices is used for determining the global solution between elements. For a onedimensional advection problem this approach is always exact as discussed in Godunov (1959) with relation to fluid dynamics equations. The same quality should hold for a 3D simulation as long as the direction vector can be precisely determined for all element faces. Conceptually, the DG method allows for discontinuities across the faces but for most subsonic CAA problems the solution should remain continuous. However, such feature is desirable for the coupled method since solution jumps may appear when receiving information from another solver. In case of strong wave reflections and scattering, the FMCAS solution is expected to be out of agreement with the near-field at least for the first update.
Every coupled simulation begins in a DG sub-domain with the following acoustic boundary conditions: (1) Full-slip condition is applied at solid walls which is deemed to be appropriate for most CAA applications. The wall boundary condition could be prescribed directly on selected cell faces by introducing a specialised treatment where the normal velocity component is removed from fluctuations. Our practice showed that such treatment is stable only for a very mildly compressible mean flow. For improved stability, the face flux can be calculated using a ghost state (Mengaldo et al., 2014).
whereV g = {û gvgŵg } T andV in = {û invinŵin } T are fluctuating velocity vectors on the ghost and inner side of the same degree of freedom respectively, andp g =p in which completes the slip wall definition for all four variables contained inq. The boundary flux is calculated via the characteristic decomposition such asH I i =H I i (q in ,q g ) at every high order point on a boundary face.
(2) The standard radiation boundary conditions are prescribed to outflow boundaries, and likewise the ghost cell method is applied to boundary faces such that the information can be received from the FM-BEM by updating the incoming flux on the ghost side. The forward coupling, DG to FM-BEM, is performed via a permeable surface placed in the region of potential flow. (See Figure 2) The surface is triangulated to meet the requirements of the FM-BEM. Then, the DG solution interpolated to the radiation surface becomes part of the BEM problem which now includes the entire scattering geometry. Thus, the FM-BEM problem can be solved for the remaining unknowns based on the partly given surface solution.
Unlike with most fine-scale turbulence resolving methods where the location of an integration surface affects the quality of CAA results, e.g. due to excessive numerical dissipation, such limitation does not apply to the wave propagation computed by DG. In general, one is free to set the appropriate scope as long as the coupling surface is in a potential region. Particular attention is given to communication between CAA regions since the ability to accurately exchange data is the key feature of coupling. The information between DG and FM-BEM can be exchanged via the fluctuating enthalpy using the Möhring-Howe acoustic analogy (Heitmann et al., 2016). Alternatively, one may choose to reformulate complex pressure and velocity values in terms of the acoustic potential based on properties of the wave equation. In the Möhring-Howe approach, the main variable of the perturbation system is the fluctuating enthalpy, given for an ideal gas:B In an isentropic region ρ = p /c 2 0 , substituting into the above equation yields, where the acoustic pressure and velocity computed in a mean flow are readily available from the DG simulation.
For an irrotational flow region, the acoustic velocity can be related to the gradient of the fluctuating enthalpy, v ≡ ∇φ = − ∇B iω (18) and upon proper gauge the above relationship infers: As discussed in Section 2.1, the Taylor phase factor is applied to the potential in the framework of the Taylor transformation,φ = Tψ . Next, the Taylor transformed acoustic potential can be derived from Equations (17)-(19) in terms of the DG variables on the coupling surface (indicated by subscript s), Also, the normal gradient ∂φ s /∂n is required. The following expression can be deduced at the coupling surface: recalling that T = e iκM 0 (x)/v ∞ = e −iκM (x) , it follows that ∂T ∂nψ s = −B s (M∇ · n)/c 0 . Then using the properties of Equations (18) and (19) in Equation (21) and substituting for the fluctuating enthalpy using Equation (17) one can obtain the surface normal derivative of the Taylor transformed acoustic potential in terms of the DG variables on the coupling surface up to linear order in M, By examining Equations (20) and (22) it is possible to see how information is exchanged by using the Möhring-Howe acoustic analogy. On one side the DG method operates withp andv variables and on the other side the FM-BEM uses the Taylor transformed acoustic potential which is derived from the fluctuating enthalpy for a surface source. The feedback information arriving to the DG boundary ghost points must contain a full velocity matrix in addition to pressure to comply with the APE system. Since velocity vectors cannot be supplied directly from the FM-BEM solution, they have to be extracted from the gradient of fluctuating enthalpy, refer to Equation (18). The gradient is evaluated as part of the FM-BEM solution. At the same time, the pressure is evaluated by means of Equation (17) using thev term that has just been derived, and thus, the mean flow term in Equation (17) can be also accounted for in the DG system.
It is important to recognise that incoming characteristics may significantly alter the DG solution and when convergence within the source region is reached for the second time, data exchange takes place once again. The iterative process repeats until the residual is sufficiently small. From experience, only one or two updates are required to arrive at an established solution which also accounts for the source correction.
In summary, one iteration cycle of the two-way data exchange depicted in Figure 2 can be described as follows.
DG FM-BEM: (1) Acoustic pressure and velocity data (p ,v ) are interpolated from a volume DG grid onto a coupling surface which consists of elements compatible with the FMCAS (FM-BEM) code.
(2) The acoustic variables on the coupling surface are cast into the Taylor transformed acoustic potentialψ s and its gradient ∂ψ s /∂n by means of Equations (20) and (22) (1) From the given FM-BEM solution forψ , the FMM plotting function (evaluation of volume data from the surface solution) is applied for rapid calculation of the potentialφ and its gradient ∇φ ≡v at the given DG outflow boundary points where Equations (18) and (19) are utilised going via the fluctuating enthalpy.
(2) The fluctuating pressure is computed from Equation (17):p = ρ 0B − ρ 0 v 0 ·v . (3) Finally, (p ,v ) data are used in the APE equations to prescribe incoming acoustic wave modes via the characteristic far-field boundary condition, and thus, resuming a volume resolving DG step.

Acoustic sources in DG
Acoustic sources can be defined on a surface or in a volume depending on the source type. For concept validation, a monopole source is placed in front of a rectangular screen. The source is enclosed by a box meshed with tetrahedra. Overall, the DG domain has five radiation faces and one solid boundary (on the face adjacent to the screen). Constant mean flow of Ma 0.2 is prescribed parallel to the solid boundary. Figure 3 depicts the acoustic potential evaluated on the BEM surface. The rectangular screen of finite thickness has a solid wall boundary condition whereas the protruded part is a permeable source surface which can be referred to as the DG / FM-BEM coupling surface. Both the screen and radiation parts of the BEM surface feature an opening but together form a common watertight surface. Any reflected acoustic waves coming from the screen do not see the source surface as an obstacle and the idea is to prescribe accurate incoming characteristics to the restricted DG domain. Notice that at first, the DG solution cannot account for acoustic shielding effects on its own since there is no information about the scattering body. The FM-BEM solution, on the other hand, can be used to predict the acoustic potential on the solid rectangular screen of finite thickness, subsequently establishing the global solution. For evaluating the surface integral in Equations (9) and (10), the FM-BEM needs pressure and velocity vectors supplied from DG which have to be recast into surface sources. More sophisticated acoustic sources other than generic monopoles can be also prescribed in the frequency version of the DG code. Channel modes can be inserted at circular boundaries closely resembling engine fan noise. This is particularly useful when a precise source definition cannot be obtained from experimental measurements, e.g. for a scaled model. The DISCO++ user parameters can take into account the definition of azimuthal and radial modes, amplitude, phase shift, inner and outer radii, local flow Mach number, frequency and propagation direction. Since noise simulations are integrated in the design study, the above source model is applied in Section 6 for simulating the propagation of fan modes inside the intake duct of a UAV model. In some cases, it may even be possible to make quick noise predictions with broadband sources, sampled on a surface in the frequency space (Reiche et al., 2017). Before proceeding with the DG / FM-BEM coupling for engineering cases the scattering of an acoustic monopole by a thin flat plate is considered.

Monopole scattering by a plate with flow
A test case performed by Lummer (2019) with wavelength λ = 0.25 and 0.5 m was revised and extended to higher Mach numbers. The selected wavelengths correspond to the Helmholtz number (He = 2π L/λ) of approximately 50 and 25 respectively with the characteristic plate length of 2 m. In the current implementation, the FM-BEM is restricted to low Mach applications only, since Ma 2 and higher order terms must be dropped in the definition of the Taylor transformed acoustic potential. The solution errors resulting from the absence of corrective terms may grow rapidly with the increasing Mach number. Thus, one objective was to verify the limit at which the solution accuracy of the FM-BEM would be no longer acceptable for noise predictions when solving engineering problems. The other objective was ensuring that the same monopole source description was used in both solvers not only in terms of spatial parameters such as position and half-width but also concerning its physical type and amplitude scaling. It is not uncommon for different methods to employ different source models so this step was essential because fine differences help to reveal and better understand problems with coupling which may not be picked up otherwise. Table 1 shows the simulation parameters. All DG simulations were executed on the DLR cluster CARA 2 using 32 nodes (4 processors, 64 threads on each node) while the FM-BEM simulations were already efficient on a local machine running just 10 threads. The number of elements and degrees of freedom of the mesh are displayed in millions. The next column contains an estimate of points per wavelength without accounting for the Doppler shift. A couple of quick calculations revealed that acoustic data were resolved on the DG mesh up to Ma ≈ 0.6 at the numerical microphones located in the close proximity of the plate. This resulted in a suitable testing range. Due to different aims and objectives to the test case performed in Lummer (2019), the FW-H surface was not used taking any uncertainty of additional data processing out of the equation. Also, in contrast to the original test case which aimed for optimal performance, this time the simulations for a different λ were performed on the same grid (depending on the solver). Hence, the grid-related parameters remained identical and only the number of iterations and the Wall Clock Time varied between λ cases due to subtle differences in solution advancing and convergence. In the FM-BEM the GMRES solver was used by default and the DG simulations were performed using the explicit RK4 scheme in order to give the reader a better feel for a typical setup and computational requirements.
A flat plate with 2 m sides and 0.04 m in thickness and the midpoint located at x, y, z (0, 0, 0) fully replicates the plate used in Lummer (2019). A compact monopole source with the half-width of 0.04 m was placed asymmetric (0.7, 0.7, 0.1) to the plate geometry close to the trailing edge of the upper surface. A constant mean flow was prescribed in the x direction. The computational DG domain was a cube with dimensions −3 < x, y, z < 3. Figure 4 shows a slice cut off-centre through the 3D mesh. The mesh was refined in the source zone and close to the surface fitting at least a couple of elements per plate thickness around the edges. Figure 5(a-c) show the contours of pressure for λ = 0.25 and the Mach numbers For a zero flow condition a perfect match is achieved for each wavelength. The agreement at Ma = 0.2 is consistent with the comparison shown in Lummer (2019) but further doubling of the Mach number resulted in a more pronounced deviation in all directions, especially for the central lobe in Figure 6(c). A further inspection of the Ma = 0.4 acoustic solution revealed that the source centre point was displaced due to the Doppler effect which possibly caused a strong trailing edge response for λ = 0.25. The effect may be negated by extending the plate in the flow direction but importantly, it shows that the FM-BEM cannot account for sources which are strongly influenced by convection. Any further increase of the mean flow speed resulted in further deviation of the FM-BEM pressures from the reference DG solution. Thus, it can be concluded that reliable acoustic predictions cannot be obtained if the flow exceeds Ma 0.4. In terms of the absolute noise level, in the region of 0.2 < Ma < 0.4, only some minor differences will be seen on the dB scale. However, it involves large risks of obtaining incorrect predictions for more complicated scattering problems, especially at high frequencies.

Coupled DG / FM-BEM
The geometry in Figure 3 was used to set up the initial coupling problem. In this problem a monopole source   with a wavelength of 0.3704 was placed at the centre point of a square plate with dimensions 2 × 2 × 0.2, precisely 0.5 m away from the plate. Figure 8(a) shows the triangulated surface where the DG / FM-BEM data were exchanged via the protruded purple interface which became part of the BEM surface. The real part of the pressure potential shown in Figure 8(b) was obtained in the outcome of the coupled simulation for a zero mean flow condition. The Re(ψ) solution appears to be fairly consistent across the interface despite the connection between surfaces being at a right angle which is unfavourable for computing the pressure gradients. Figure 8(c) shows the total pressure plotted for the cross-section which gives the reader a view of the scattered field. Notice a few hot   spots which appear at the scattering surface. This side effect is caused by the FMM plotting having to evaluate the solution very close to the surface inside the near field of some fast multipole cube. Initially, by solving the problem without a background flow allowed checking the scattering pattern for symmetry. See Figure 9 which shows the SPL contour comparison plotted on a large (20 × 20) plane six screen lengths below the source. At that point the coupled solution underwent three cycles of back and forth data exchange and concluded to be fully converged. The reference FM-BEM solution shown on the right consisted of a point monopole source distributed over a sphere with a radius corresponding to a half-width of the DG source. One of the objectives was to ensure that the coupled approach can provide accurate noise shielding predictions, on the same level with the DG / FM-BEM comparison presented in Figures 6(a) and 7(a), for the proposed method to be taken a step further.
The DG simulations were performed up to numerical microphones similar to monopole scattering simulations presented in the previous section. The p RMS directivity data shown in Figures 10(a) through 10(c) were sampled at a radius of 2 m from the centrepoint of the plate using 360 numerical probes. The probes were rotated around the vertical axis and aligned with the flow direction parallel to the screen surface. Figure 10(a) shows that all three simulations are in perfect agreement when computed without potential. In Figure 10(b) for the Mach number, Ma = 0.2, the most pronounced difference between solutions is visible along the 180 • axis where the acoustic pressure amplitude of the FM-BEM simulation is significantly under-predicted in comparison to the reference DG and coupled solutions which both came to a good agreement. The FM-BEM simulation was repeated with an ideal point source which resulted in exactly the same behaviour. After revisiting equations for the Taylor transformed potential, (see Equation (23) for the acoustic pressure which is at the core of the transformation), it follows that the mean gradient term, ∇ , which is the local Mach number vector becomes zero when evaluated along the source axis normal to flow. For compact sources, the correction along 180 • is simply absent meaning that the amplitude and shape should correspond to the case without potential. The central lobe obtained for the FM-BEM in Figure 10 shadow side (at 0 • ) which is not visible here because all noise was shielded.
(23) For the coupled case, the source surface was about 3 wavelengths wide which led to the correct contribution at 180 • at a far-field observer using the same formulation. The problem was not encountered in the previous test case because a monopole source was placed at the edge of the plate, and thus, always positioned at an angle relative to the far-field microphones. Notice the upstream radiation from a monopole source being visible in all simulations with added potential where the flow vectors are aligned with the positive y-axis (pointing at 90 • ). Figure 10(c) continues the trend for Ma = 0.4. Here, the applied Taylor transformation under-predicted the upstream radiation by a big margin primarily due to divergence from its validity bounds which is consistent with results obtained for the first validation problem. In fact, when the Mach number correction was disabled, the upstream lobes (pointing towards 240 • ) of the coupled problem overshot the directivity obtained by DG for Ma = 0.4. This confirms that the transformation with added potential has a significant effect on the radiation amplitude.

Coupled simulations for a UAV intake noise predictions
The original UAV configuration was developed in the framework of NATO STO Task Group AVT-251 a couple of years ago. Their goal was to perform aerodynamic optimisation of the UAV model through the application of CFD methods. The AVT-318 task group worked independently on low noise aero-acoustic turbofan design based on the previous research started by the group AVT-233 whose mission was to explore the aero-acoustics of engine installation on military air vehicles. The collective effort led to the design and construction of the current UAV MULDICON model. The DLR-F24 MULDICON features a new design for intake and exhaust ducts but the original planform proposed by the AVT-251 which was later modified by FOI (Swedish Defence Research Agency) remained untouched. Figure 11 presents an overview of the model focussing on the intake and exhaust ducts. The T-junction was designed to provide a controlled high-velocity intake flow in the low-speed acoustic (DNW-NWB Braunschweig) wind tunnel.

RANS simulation
In the initial phase of the CAA study discussed in this paper, the focus was solely on the intake radiation noise. Figure 12 shows the contour slice of the Mach number for the DLR-F24 MULDICON with a realistic jet engine inlet profile. The RANS simulation was provided by Airbus Defence and Space (ADS) and the free flight conditions were chosen close to those which can be replicated for the wind tunnel tests. From the noise modelling view point this also meant that the exterior flow of Ma = 0.175 is well within the scope of the FM-BEM tool. It can be seen that the Mach number quickly rises as the flow enters the intake duct reaching a peak of 0.75 at the inlet of a virtual fan.

Intake radiation
In the acoustic simulation, the initial step consisted of generation and propagation of channel modes through the duct using the DG method with and without the background RANS flow solution. A spherical coupling surface was placed over the inlet for restricting the DG volume domain to the region of a highly non-uniform flow. Figure 13 shows the triangulated FM-BEM geometry with the attached coupling surface. According to the theory of BEM, radiation surfaces can be immersed in potential flow at most. However, this condition cannot be satisfied for a small portion of the coupling surface where it must naturally penetrate a boundary layer to form a closed surface with the UAV aircraft exterior. Also, since the FMM plotting must be used in the two-way coupled mode for performance, it can possibly cause accuracy problems at any DG point located close to the scattering body, i.e. in its 'near-field'. The second problem can be eliminated via careful integration treatment and the significance of the first problem can be assessed by comparing feedback data for the same simulation performed without the flow.
In the wind tunnel experiments, a laser pulse will be used to excite rotating modes inside the intake duct at the fan position. The numerical results for two frequencies of interest, 5 and 10 kHz, are provided as a proof of concept. Figure 14(a,b) show the real part of the total pressure fluctuations computed inside the duct for the 5kHz rotating modes. It is clearly visible that the radiation is heavily suppressed by a strong flow inside the duct to an extent where the external contours can no longer fit well on the same scale. Also, notice that the noise pattern is altered. It was ensured that the grid resolution was not a limiting factor through performing a convergence study, especially for the cases with flow. One has to bear in mind that naturally, some modes cannot propagate due to cut-on / cut-off condition imposed by the flow. It will be part of the CAA design study to obtain reliable predictions for complicated shape ducts with a varying cross-section before conducting an experiment. Figure 15(a,b) show the comparison at 10 kHz. The trend is similar for both frequencies where a much stronger radiation is found without the added flow. For the current UAV design, the noise radiates predominantly upstream in the direction towards the ground which is not optimal and highlights the importance of aero-acoustic design. The amplitude of the signal tends to decrease with increasing frequency. Also,    one can expect forward radiating fan noise to decrease with the increasing Mach number where at cruise conditions, Ma = 0.8 and at high altitudes the noise signature is expected to be negligible. Therefore, for the future acoustic design optimisation of the UAV configuration the focus remains on improving our noise prediction methods at low Mach numbers. This corresponds to takeoff and approach conditions where noise can be detected on the ground. The numerical coupling is formulated for the noise sources in a background flow typically encountered during this critical part of the mission. Figure 16 shows the real part of the total pressure for the coupled FM-BEM surface at 10kHz without the added mean flow and the noise footprint plotted on a large plane which represents the ground. The UAV model was rotated to get a better view of the coupling interface. A smooth transition of the surface pressure from the radiation to the scattering surface is a good indication that two methods communicated the data correctly. Due to a highly directed pattern of a forward radiating intake noise, the two-way coupling did not result in any significant correction to the DG solution. Yet big savings were achieved with the coupled method in comparison to potentially computing the DG propagation in a domain which has to include the UAV if not distant observers. Arguably, using a FW-H surface in place of the coupling surface would have resulted in similar noise predictions for this configuration. However, the idea of applying global boundary conditions on a highly restricted DG domain via performing the two-way coupling must be superior for computing accurate noise scattering and shielding by large geometries.

Conclusion
A medium fidelity technique that is based on strong coupling of two advanced CAA methods has been presented to accurately predict sound shielding for large Unmanned Aerial Vehicles (UAV). First step consists of the acoustic wave propagation in a highly non-uniform flow using a volume resolving method. Then the radiated waves are sampled on a coupling surface. The acoustic scattering produced by the geometry of the UAV is accounted for by the FM-BEM that has been proven to be very efficient for determining surface solutions. Depending on the configuration, performing a one way coupling may be insufficient for accurate noise predictions due to absence of feedback provided to the sound source region. Thus, the backward data exchange takes place on the DG outflow boundaries which is equivalent to prescribing global boundary conditions to the entire DG domain. The information is exchanged via a fluctuating enthalpy and the incoming characteristics introduce a change into the volume region. A background flow available for the volume region can be based on a RANS solution which naturally accounts for boundary layers.
In the current framework, all CAA simulations were performed in the frequency space where a new converged state was quickly reached from a 'virtual re-start' where feedback data were prescribed on the far-field boundary. The cycle was followed by another DG / FM-BEM exchange. It has been determined that a fully converged solution results after just a couple of data exchanges even for cases with strong wave reflections of a scattering geometry.
The noise prediction technique described in this work has been developed with the objective to quickly evaluate new design concepts. The future work would involve application of the coupled method for simulating fantone shielding for a full-scale supersonic aircraft. It can be challenging to optimise a supersonic nacelle for takeoff and approach subsonic phase. Currently, the method is limited to predictions in a free-stream flow of up to M ≈ 0.3 in the FM-BEM domain. This limitation comes from a low Mach number assumption used in the formulation of the BEM equations. It appears that there are no simple ways to overcome this restriction, also, bearing in mind that the FMM feature is very desirable and has to be retained for efficacy when solving engineering problems.