On the prediction of noise generated by urban air mobility (UAM) vehicles. II. Implementation of the Farassat F1A formulation into a modern surface-vorticity panel solver

For more than a decade, the emergence of urban air mobility (UAM) and distributed electric propulsion (DEP) vehicles^1–3 including unmanned aerial vehicles (UAVs) has given rise to a new platform to consider in the aerospace industry worldwide. New designs have surfaced that augment lifting capabilities using a large number of electric-motor-driven propellers. Naturally, modern DEP vehicle concepts aim to realize favorable interactions between the airframe aerodynamics and the multiple electric propellers that are strategically integrated into the airframe to achieve hitherto unrealizable gains in performance and energy consumption. Due to its novelty, the underlying paradigm shift in vehicle characteristics warrants new demands for an array of predictive capabilities. More specific to the present work, the intended urban environment for this class of vehicles places strict limitations on acoustic radiation. The corresponding limitations must be addressed at least preliminarily in the conceptual phases of design to ensure the viability of the final product. With the physics of incorporating aeroacoustic predictive tools into a robust flow solver being the main objective here, a companion article addresses the physical definitions and interpretations of the acoustic metrics that can be evaluated to support the UAM design cycle.⁴ 

Due to the appreciable growth in civilian populations and the number of airports and aircraft globally, the ability to predict and mitigate aircraft noise stands as a problem of tremendous interest. It continues to occupy center stage in the aircraft mobility industry in general, and the UAM research community in particular. Studies focused on the latter include, for example, those by Jia and Lee,⁵ Lee and Lee,⁶ Gwak et al.,⁷ Smith et al.,⁸ Krishnamurthy et al.,⁹ Jeong et al.,¹⁰ Dbouk and Drikakis,¹¹ Ko et al.,¹² and Wang et al.¹³ To engage in the prediction of noise from a given UAM vehicle, one must begin by seeking to understand the relative contributions of noise-generating sources. To this end, and in the spirit of illustrating the wide range of sources that must be accounted for, particularly in what concerns the UAM vehicle category, it may be helpful to refer to the classification of noise sources according to Li and Lee,¹⁴ which is outlined concisely in Fig. 1.¹⁵ 

It is well accepted that the most significant sources of noise for propellers with subsonic tip speeds consist of the tonal noise components of loading and thickness types.^16–18 These two primary sources of noise will therefore constitute the main focus of the present investigation and integration effort. In this process, a methodology will be established that can be readily extended to perform multi-rotor DEP concept simulations to the extent of producing a conceptual-phase aeroacoustic toolbox within the solver. From a practical perspective, this choice of formulation will be well suited to handle the conceptual design problem as it strikes a reasonable balance between modeling fidelity and computational tractability, especially that it can be embodied into a dependable and fast unsteady solver.^19,20 The latter consists of a novel surface-vorticity unstructured solver that is developed in-house. It will be relied upon to generate the unsteady aero-propulsive loads on the UAM or DEP vehicle in a flight mode with spatially and temporally varying velocity fields and propeller operating conditions.

To address the overarching requirement for efficient prediction of UAM noise, the surface-vorticity unstructured panel method used by Ahuja and Hartfield¹⁹ and DiMaggio et al.²⁰ will be coupled with a simplified acoustic formulation based on the Farassat 1A solution of the Ffowcs Williams–Hawkings (FW–H) equation.²¹ The F1A formulation represents a viable alternative solution to the FW–H equation for thickness and loading noise computations by surface integration through the vorticity-flow solver.²² The utilization of surface-vorticity and vortex methods in this context offers unique advantages; this has been demonstrated on several occasions by Maskew,²³ Dvorak et al.,²⁴ Maskew,²⁵ Katz and Plotkin,²⁶ Wu et al.,²⁷ Ahuja et al.,²⁸ Ahuja and Hartfield,¹⁹ Faure and Leogrande,²⁹ Faure et al.,³⁰ and Fenyvesi et al.³¹ As for the solver itself, several validation studies have been undertaken at the National Aeronautics and Space Administration (NASA) Langley Research Center³² as well as the UAM industry.^19,20,33

As part of prior research activities, the flow solver has been developed to make use of surface vorticity on an unstructured surface mesh and predict loads with attached flow.^19,20 Within its framework, aerodynamic loads are computed by shedding vorticity from the object to be analyzed and then tracked using the fast multipole method (FMM). The general implementation of FMM in any fluid-flow problem revolves around coalescing far-field or potential inviscid sources and doublets into a single body for evaluation of the net effect on a near-field surface.¹⁹ The algorithm evaluates the inductive effect of near-field bodies identically to the $O(n2)$ implementation although it converts the large number of far-field bodies into a handful of clustered objects by taking into effect the sum total of their field strengths. Within the solver, FMM implementation condenses far-field vorticity into localized three-dimensional point doublets encompassing the net vorticity strengths of all of the coalesced far-field vorticity panels. Once the FMM near-field and far-field thresholds have been established, a spatial tree is constructed around the three-dimensional geometry to allow for the demarcation of these thresholds and to identify the near-field neighbors.

At the root of the present solver, a doublet-based vorticity scheme has been developed and implemented into the solver.¹⁹ In this context, the spherical harmonics of degree n and order m of the underpinning formulation can be expressed as follows:

where θ and $ϕ$ denote the azimuthal and polar coordinates of a spherical reference frame. In the above, $Pnm(x)$ represents the Legendre functions given by

and

Let us now consider a three-dimensional point doublet with strength μ. As described by Katz and Plotkin,²⁶ it is well known that the potential induced by a point doublet in a generalized coordinate system and its relation to a point source of strength σ can be expressed using the following expression:

where $n̂$ stands for the outward pointing surface-normal unit vector comprising $n̂i$ components and r refers to the radial distance from the source. To make further headway, the numerical solver follows the same mathematical strategy developed by Cheng et al.,³⁴ which is outlined here for clarity. Suppose that N charges of strengths ${q1, q2, …, qN}$ are located at points ${X1$⁠, X₂, $…, XN}$⁠, with spherical coordinates (⁠$ρ1, α1, β1$⁠), respectively. Suppose further that the points ${X1$⁠, X₂, …, $XN}$ are located inside a sphere of radius a that is centered at the origin. Then, for any point $X=(r,θ,ϕ)$ with r > a, the potential $Φ(X)$ generated by the doublet strengths ${q1$⁠, q₂, …, $qN}$ may be written as follows:

Here $Mnm$ stands for the multipole term that combines the far-field potential functions of the individual doublet singularities whose spherical harmonics of degree n and order m is given by $Yn−m$⁠. As for the radial distance, r, it extends from the center of the multipole sphere enclosing the doublet singularities to the far-field point X where the potential $Φ$ is being evaluated.

In the context of this study, both laminar and turbulent boundary-layer models are judiciously implemented. A boundary-layer transition model is also incorporated. All of these models translate two-dimensional formulations along the on-body surface streamlines. Furthermore, due to their computational efficiency, these integral methods can be applied to any general three-dimensional wall boundary except in the regions involving crossflows. The laminar boundary-layer technique used in the solver corresponds to the standard two-parameter model of the Thwaites integral method paired with the momentum integral equation.^35,36 It is given by

where U, θ, x, and ν denote the freestream velocity, momentum thickness, longitudinal streamline coordinate parallel to the surface, and kinematic viscosity, respectively. An integral boundary-layer model for compressible turbulent motion is also embodied within the inviscid-flow solver.¹⁵ The final turbulent boundary-layer model incorporates improvements and optimizations over that of the original model developed by Standen;³⁷ nonetheless, it conceptually follows the formulation by Standen quite closely, thus leading to effective applications to subsonic, turbulent, and compressible motions along the on-body streamlines. At its heart, the method remains essentially two-dimensional and similar to the laminar model described previously, albeit extendable in a quasi-three-dimensional manner inside the solver. Moreover, in evaluating fluid properties outside the boundary-layer region, the flow may be taken to be isentropic and compressible within the subsonic regime. As such, the primary differential equations are those associated with the turbulent compressible momentum thickness, $θc$⁠, as well as the compressible shape factor, H_c. These follow from the coupled equations:

and

where

In the above equations, C_f , θ_i , H_i , H_TR , $Hδ−δ*, Fδ−δ*$⁠, and $Mae$ represent, sequentially, the skin friction coefficient, the incompressible momentum thickness, the incompressible shape factor, the transformed shape factor, the entrained shape factor, the nondimensional mass entrainment rate, and the Mach number at the boundary-layer edge. Moreover, T₀ , T_aw , T_w , T_R , and T_e refer to the total stagnation temperature in the far field, the adiabatic wall temperature, the wall temperature, the reference temperature, and the static temperature at the boundary-layer edge, consecutively. As for δ and $δ*$⁠, they stand for the compressible boundary-layer disturbance and displacement thicknesses, whereas μ₀ and μ_R denote the dynamic viscosity at T₀ and T_R, respectively. Given the number of unknowns, three auxiliary relations are still required to achieve closure. These include the definition of the Mach number, the derivative of the Mach number along a streamline, and the definition of the shape factor. As shown by Standen,³⁷ these can be retrieved from the following:

and

The additional parameters, U_e , c_e , γ , g₀ , R₀ , and ρ₀ allude to the far-field velocity at the edge of the boundary layer, the sound speed at the edge of the boundary layer, the ratio of specific heats, the gravitational constant, the individual gas constant, and the density at T₀, respectively. In practice, the resulting equations can be readily solved numerically for all relevant boundary-layer characteristics using, for instance, a conventional fourth-order Runge–Kutta scheme along the on-body streamlines.

It can thus be seen that with the integration of modern viscous models, several computational enhancements can be incorporated into the solver. This can be accomplished while still maintaining an $O(n log n)$ scalability to the performance of the boundary-layer algorithms. To this end, spatial octree algorithms that mimic the existing FMM implementation within the inviscid solver can be extended to the viscous near-wall regions. This extension leads to rapid near-field sorting of points that are closest to or inside a viscous-flow region. The corresponding octree-based volume grid, termed the viscous spatial tree (VST), can be further refined to capture the three-dimensional volume of the boundary-layer region extending above the surface of the geometry. Following this corrective implementation, the VST increases the memory footprint of the otherwise inviscid solver by less than 5% for all cases considered. However, the performance benefits are substantial: in lieu of an anticipated $O(n×m)$ algorithm, where m denotes a very large integer corresponding to the total number of surface mesh faces, one essentially arrives at an appreciable reduction to $O(n log n)$⁠.¹⁵ 

Apart from these boundary-layer refinements, a robust model for computing attachment lines and critical points on the inviscid surface velocity fields is integrated into the solver. This particular step enables us to systematically retain the flow coupling with the viscous boundary-layer models described previously. In practice, the laminar and turbulent boundary-layer formulations are embedded within the inviscid solver via displacement of the inviscid boundary by an equal amount to the displacement thickness of the local boundary layer. The latter indicates the extent to which the surface has to be displaced to permit the same flow rate as the viscous motion while using an inviscid velocity profile. The conventional displacement thickness concept can, therefore, be used to predict the required displacement of the inviscid boundary. In this work, the inviscid boundary displacement is simulated by a transpiration flow boundary condition that is added to the inviscid Neumann-type constraint normal to the mesh face. The magnitude of the transpiration velocity, V_N, can be computed rather straightforwardly using the momentum flux equation,¹⁵ which is given by

where $δ*$ denotes, as usual, the displacement thickness based on the integral boundary-layer computations; on the other hand, U_e stands for the local streamline velocity directly outside the boundary layer, and s represents the direction along the local surface. Equation (14) can be computed for the transpiration velocity V_N along the on-body streamlines for each mesh face using the standard finite-difference numerical marching schemes.

Over the years, NASA has created, or has had access to, a variety of high-order acoustics formulations for use with both potential-theory solvers and conventional Navier–Stokes CFD solvers. Tools such as WOPWOP, WOPWOP+,^38,39 ANOPP2,⁴⁰ ASSPIN,⁴¹ and others^42,43 have provided very powerful acoustic analysis capabilities. In most cases, however, they have lacked a user-friendly interface or a convenient post-processing toolbox to encourage their adoption by the nonspecialists. Gradually, these tools have been refined to incorporate the high-order nonlinear terms into the acoustics framework. To clarify, it may be helpful to illustrate through superposition the various acoustic components that prescribe the total acoustic signal at an observer location at a given time. One can calculate

where $p′$ represents the acoustic pressure perturbation, x denotes the observer position relative to the source, and t refers to the observer's time. In the above equation, the term $p′Q(x,t)$ relates to the acoustic quadrupoles, while the $p′NL(x,t)$ term captures the correction due to nonlinear effects. The first two terms, $p′T(x,t)$ and $p′L(x,t)$⁠, correspond to the thickness and loading components of the acoustic signature at the observer's location, respectively.

It should be noted that the noise due to compressibility effects is typically accounted for in a third term, namely, the quadrupole noise source. In the present implementation, the quadrupole noise source is not considered because it reflects a volume source that requires a grid-based Navier–Stokes flow solver; as such, it will be quite expensive to evaluate in early design workflows. In fact, the volume integral of the quadrupole sources that arises in the nonlinear region outside of an acoustic control surface presents a major challenge and adds significant complexity to the acoustics models. However, if one ignores the quadrupole and nonlinear terms, the resulting framework reduces the necessary calculations by one order of magnitude, specifically, by lowering the order of the acoustic solution. Similarly, if the nonlinearities on the control surface are ignored, the acoustics computations assume a solution to the linear wave equation on the surface. From a practical perspective, and for both DEP and UAM flight vehicles, the nonlinearities can be reasonably ignored, especially in early design stages. At the outset, the total signature equation turns into the sum of the thickness and loading terms,

On the one hand, the thickness noise source $p′T$ accounts for noise due to the displacement of the fluid by the finite thickness of the body. On the other hand, the loading noise source $p′L$ accounts for noise due to loading and change of loading on the body. The typical starting point for any numerical noise prediction model attempting to solve for the thickness and loading terms is one of the various forms of the solution to the Ffowcs Williams–Hawkings (FW–H) equation.²¹ In the present work, the FW–H equation for a permeable surface is employed. This is known as Formulation 1A and yields the following equations for the thickness and loading terms:²² 

and

where

In this formulation, the subscript i in tensor notation implies an inner product between two vector components. Other terms such as c, $dS$⁠, v_i, and $ℓi$ allude to the average speed of sound, the infinitesimal element of surface area over which integration is carried out, the local fluid velocity component, and the local force intensity component acting on the fluid, specifically $Pijn̂j$⁠, where P_ij denotes the compressive stress tensor.

In the present solver, these equations are integrated using the time-dependent aerodynamic data, and the corresponding derivatives are generated by the unsteady solver. The rotor blades are discretized into panel elements within the solver, each of which being an acoustic source. On the one hand, the loading noise is calculated based on the loading and velocity of each element. On the other hand, the thickness noise term is retrieved from v_n, the blade's local normal velocity. The loading and thickness noise contributions of all elements are then collected to deduce the total acoustic noise signal at user-defined measurement points.^17,18

It may be instructive to note that, in the loading and thickness expressions, integrands with $r−1$ dependence denote far-field terms that decay rather slowly, whereas those with $r−2$ dependence constitute near-field terms that dissipate quite precipitously. Although the integrands in Formulation 1A are somewhat more complicated than in Formulation 1, a numerical differentiation of an integral is no longer required. This is especially beneficial for cases where the observer is moving.

To facilitate the implementation of Formulation 1A, the retarded time computation scheme chosen for this work will be based on the source-time-dominant model described, for example, by Brentner and Farassat.¹⁸ This model considers the source time as the primary or dominant time. Rather than select the observer time in advance, one can choose the source time for a panel i at each panel center and determine when the signal will reach the observer. If the observer x happens to be stationary, then one must have

In the above equation, r_i denotes the radial distance from the acoustic source (located at $yi$⁠) to the observer (at x), whereas t and τ refer to times measured at the observer and source locations, respectively. In most cases, this time projection can be computed immediately. Otherwise, one can extract the root of

The determination of t even for the latter case is easier than finding the retarded time because the observer motion is usually simple to track; in this event, the solution for t may be found analytically rather than by iteration. A sequence of source times (i.e., the times at which the source strength is available) ultimately leads to a sequence of unequally spaced observer times. This panel time history can then be interpolated to provide the proper contribution of each panel to the net acoustic signal as seen at the observer location. Interpolation in time is subsequently undertaken to consolidate the signal contributions originating from all source panels at the same observer times.

This foregoing algorithm has the unique advantage that a retarded-time calculation is not necessary per se and that the discrete time-dependent input datasets do not need to be interpolated. This characteristic can be very beneficial, especially when the surface-vorticity solver itself is being relied upon to generate the input data internally and seamlessly. Another computational advantage of the source-time-dominant algorithm can be ascribed to the solution process being inherently parallel; in fact, the algorithm proves to be an excellent candidate for parallel computing. Performance tests performed by Metzger³⁸ confirm that the source-time-dominant algorithm requires significantly fewer operations for a maneuvering vehicle prediction. Finally, although an effort has been made to investigate quadrupole term implementation, it has not been further pursued based on the knowledge that most UAM vehicles are not expected to fly in the regimes where quadrupole terms become appreciable or even measurable.

Acoustic volume sections are planar grids of stationary microphones placed in an arbitrary orientation, size, and resolution around a vehicle being simulated using an unsteady surface-vorticity flow solver. The volume section can be used to generate a contour of sound pressure–time history that can then be exported and analyzed for acoustic volume field effects rather than point effects for a single stationary observer. To proceed, a representative aircraft is selected, particularly, the NASA Revolutionary Vertical Lift Technology (RVLT) Tiltwing electric Vertical Take-Off and Landing (eVTOL) vehicle. The corresponding geometry is based on the Open Vehicle Sketch Pad (OpenVSP) file of the RVLT concept. As shown in Fig. 2, this particular RVLT vehicle is powered by eight 5-bladed propellers operating at 1420 RPM and a tip Mach number of 0.5. The vehicle is simulated in the forward-flight mode at a zero angle of attack and a cruise speed of 20 m/s. In Fig. 2(a), 24 536 surface panels are shown around the vehicle. In practice, however, one may prescribe symmetry about the midplane of the aircraft and reduce the number of simulated panels by half. Although different aerodynamic structures and loads can be deduced from the present surface-vorticity solver,¹⁹ an isometric view depicting vorticity streaklines is provided in Fig. 2(b). The main features that can be inferred therein consist of the induced velocities from the propellers in the forward-flight mode, which are captured using the unsteady, time-accurate flow solver. The interaction of the wake from the forward propeller with the rear propellers is also visible along with the minor distortion of the vortex field. The propellers are modeled in this simulation using thin-blade surface representations that are consistent with vorticity theory.²⁷ The actual simulations are run for a total of 144 iterations in time, with an increment of 0.002 s, thus leading to a total simulation time of 0.288 s. One of the main advantages of using this solver is that, computationally, only 500 internal iterations are typically sufficient to achieve convergence while using modest resources that include, in the present case, five physical cores and a standard workstation.

Forthwith, examples of acoustic volume sections surrounding the RVLT are produced and provided in Fig. 3. Therein, the outer radius of the acoustic volume section is set at 12 m, and the grid resolution is taken to be 40 × 40 in the radial and azimuthal directions, respectively.

Figure 3 clearly illustrates the nature of the acoustic pressure perturbation field around the vehicle using an orthogonal XYZ coordinate system; within this reference frame, the horizontal X and Y axes run parallel to the fuselage and wings, respectively; here, Z stands vertically, i.e., normally to the planform area. The acoustic pressure $p′$ values are based on the oscillatory component of the signal, which is computed through Eq. (16). In the XY plane [top view in Fig. 3(a)], the strong signal characteristics are visible in the section located 3 m below the vehicle centroid. By the time that the signal reaches 20 m below the vehicle, the signal becomes quite faint [Fig. 3(b)]. Also visible are the strong signals engendered near the tips of the propeller rotor disks, as seen in the YZ and XZ section fields of Figs. 3(c) and 3(d). Note the phase change of the signal strength upstream and downstream of the rotor disk. Interestingly, an acoustic noise “bubble” develops around the vehicle near field in the YZ plane, namely, 5 m upstream of the vehicle centroid, where it just touches the nose of the vehicle [Fig. 3(c)]. These visual characteristics are complemented by the side view of the aircraft, which is taken 5 m to the left of the vehicle [Fig. 3(d)]. At this juncture, one may recall from Sec. II A that the acoustic shielding effects off the fuselage are not presently modeled. This explains, in part, the clean nature of the noise bubble around the vehicle. In fact, acoustic scattering constitutes another mechanism that is not implicitly captured by the present implementation. In the companion study,⁴ the fast Fourier transform (FFT) functions for the present acoustic fields are discussed along with the metrics for determining the overall sound pressure level (OASPL), sound exposure level (SEL), power spectral density (PSD) spectrum, proportional frequency-band spectrum (PBS), perceived noise level (PNL), and effective perceived noise level (EPNL).

Acoustic spheres surrounding an object to be modeled constitute a direct extension of the acoustic volume section concept, although they seem to allow for a more detailed volumetric signal field analysis. The acoustic sphere is conceived as yet another user-controlled step for post-processing in the surface-vorticity solver. As in the case of the acoustic volume sections, the acoustic sphere can be placed in any orientation, size, and resolution around a vehicle being simulated using the unsteady flow solver. Moreover, as done previously with volume sections, a sphere can be used to generate vivid contours of sound pressure–time history that can be subsequently exported and analyzed for the acoustic volume field effects rather than just the point effects for a single stationary observer.

Figure 4 provides two examples of the acoustic spheres generated around the representative RVLT vehicle in the forward-flight mode. In this example, all eight 5-bladed propellers are in synchronous operation, and the images provided display two snapshots taken at two particular instants of time, t = 0.04 and 0.2 s. Acoustic pressure values (in Pa) may be inferred from the surface of the sphere. Signal data can, in turn, be readily computed at the centroid of each spherical surface panel as a function of time. This enables us to project the complex vehicle signal characteristics onto a simplified spherical surface, which can then be post-processed rather straightforwardly for far-field signal interpretation. Note that, through Fig. 4-type graphical representations, one can achieve a very dynamic visualization of the change in acoustic pressure as a function of time on the surface of the sphere, particularly, as the tips of the propeller blades make their closest passes to different segments of the sphere at various observer times. Additional examples of the acoustic spheres are provided in Sec. III.

Since one of the main objectives of this study is to demonstrate that the methods chosen for implementation are reasonably accurate and viable for use in an early design analysis, it is useful to validate the numerical implementation through comparisons with recently acquired experimental measurements, such as those obtained by Litherland et al.;⁴⁴ these correspond to the high-lift propeller (HLP) used during takeoff and landing of the NASA X-57 Maxwell DEP aircraft.

The NASA HLP validation case involves the use of acoustic data collected at the NASA Langley Research Center for the X-57 test bed.⁴⁴ The baseline aircraft for the X-57 consists of a two-engine, gasoline-powered Tecnam P2006T. The modified aircraft is equipped with a smaller wing that is augmented for takeoff by a row of folding propellers along the leading edge, namely, to help induce the blown wing effect. The corresponding folding propellers, which are tested by Litherland et al.,⁴⁴ will be modeled here using the surface-vorticity and acoustics solvers described in Secs. I and II A, respectively.¹⁹ 

Although the X-57 HLP is described in detail by Litherland et al.,⁴⁴ it may be summarized as a five-blade propeller that is designed to produce a relatively uniform-induced velocity aft of the propeller according to the method used by Patterson et al.⁴⁵ The blades employ a constant MH114 airfoil, which is modified to include a 0.020-in. blunt trailing edge thickness by relofting its upper surface to the extent of promoting ease of manufacturing. The baseline operating conditions and dimensions of this particular HLP are provided in Table I. The actual design flight conditions that are tested include three particular flight velocities: 29.8, 38.6, and 46.3 m/s. The HLP blades are manufactured using injection molding and a long carbon fiber-reinforced polyphthalamide (PPA) thermoplastic composite.⁴⁴ 

Both the right and left turning propellers, which are dubbed right hand (RH) and left hand (LH), are tested in the Langley Low-Speed Aeroacoustics Wind Tunnel. To accomplish this, the individual propellers are investigated separately, with acoustic data being recorded using a ceiling-mounted microphone array.⁴⁴ In this process, the blade passage frequency (BPF) sound pressure levels (SPLs) are recorded as a function of the observer angle. Unexpected discrepancies between the experimentally measured acoustic signals of the left- and right-turning propellers are discussed in some detail by Litherland et al.⁴⁴ 

To produce comparative data, simulations of the X-57 HLP are performed using the methodology outlined in Sec. I. Using the present surface-vorticity solver, the velocity contours on the HLP are computed in the high-performance takeoff condition using 4800 RPM and a freestream velocity of U = 38.58 m/s. The resulting contour lines are illustrated in Fig. 5. Subsequently, two aero-propulsive performance metrics of the HLP, the thrust force and torque coefficients, are computed and compared with NASA's wind tunnel experimental data and XROTOR simulations, which are reported by Litherland et al.⁴⁴ These two metrics are carried out for U = 29.8, 38.6, and 46.3 m/s and shown in Figs. 6 and 7, respectively.

First, in Fig. 6, the thrust force coefficient, C_T is plotted as a function of the propeller advance ratio, J. As one may infer from the graph, excellent correlations can be observed between the data generated using the present approach and the XROTOR simulations at relatively low advance ratios (J <1).⁴⁴ Although both numerical approaches slightly overpredict the actual measurements, the discrepancies seem to fall well within the experimental error margin of ±5%. Moreover, a closer agreement with the LH turning experimental data (hollow symbols) may be seen across a wide range of advance ratios. At higher advance ratios (J >1), however, the present solver seems to further overshoot the thrust force coefficient obtained experimentally, and this behavior may be attributed to the solver's underprediction of flow separation effects.

Along similar lines, the torque coefficient may be calculated and plotted in Fig. 7 as a function of J. Clearly, the predicted torque stands in excellent agreement with XROTOR simulations. It also agrees quite favorably with the experimentally obtained measurements, especially those corresponding to LH turning propeller operation (hollow symbols). The agreement with the experimental data slightly deteriorates for J < 0.8, where the predicted values begin to exceed those obtained experimentally, including both LH and RH cases.

Having validated the aero-propulsive predictive capability of the present solver, attention is now turned to its aeroacoustic framework, which is based on the approach outlined in Sec. II A. It is further described in the companion study by Little et al.⁴ To this end, the present solver is used to generate acoustic pressure–time histories at points that spatially coincide with the microphone array locations identified by Litherland et al.⁴⁴ Acoustic predictions generated by the numerical solver are then processed using an assortment of acoustic metrics.⁴ This effort begins by centering the pressure time signals about their mean values to accurately prescribe the corresponding pressure perturbation. This may be accomplished using a recursive subtraction of the mean pressure values until the time-averaged pressure is reduced to a practically insignificant value by insisting on

The above expression ensures that time-independent deviations in the pressure, which may be present in the original signal, do not contaminate the acoustic spectral data. In this manner, once the pressure preprocessing is completed, the spectral density function $pf2(f)$ may be computed using a discrete Fourier transform (DFT), viz.,

where p_n represents the Fourier series coefficient associated with each frequency f_n. The series coefficients may be approximated by scaling the coefficients of the discrete Fourier transform using

where $p̂n$ denotes the nth discrete Fourier transform coefficient. As for the SPL, it may be calculated using the standard expression:

where L_p denotes the SPL, $(p2)av$ stands for the average squared pressure value in a specific frequency band, and $pref$ returns the universal reference pressure of $2×10−5$ Pa. As usual, the mean squared pressure $(p2)av$ may be calculated by integrating the spectral density function:

where $(pb2)av$ yields the average squared pressure value over a frequency band that remains bracketed by the lower and upper limits, f₁ and f₂. The overall SPL may be calculated from a spectral density integral over the entire frequency range, while band-limited SPLs are computed over smaller frequency ranges.

It should be noted that the resolution of the spectral density function remains limited directly by the DFT band thickness, $Δf=1/T$⁠, where T refers to the period of the total signal. Due to this constraint, the levels for individual frequency bands can be approximated using

where $pf2(f)$ denotes the value of the power spectral density band containing the frequency f. Similarly, $(Δf)ref$ refers to the bandwidth of the power spectral density, previously expressed as $Δf=1/T$⁠.

In what follows, the results of the present acoustic analysis as well as those generated by NASA's ANOPP2 computations and experimental measurements are summarized in Fig. 8. This is accomplished by evaluating the BPF SPL as a function of the observer angle in the $30°–150°$ range. This range follows the same geometric specifications used in the NASA experiments.⁴⁴ Thus, in conformance with the latter, measurements of the acoustic signals are showcased for both RH (hollow squares) and LH propellers (solid squares) alongside the NASA ANOPP2 simulations (broken line).⁴⁰ The results based on the present solver are overlaid using a solid line.

As one can readily infer from Fig. 8 that the low-frequency tonal noise signals are accurately captured by the present simulations across most of the observer angle range; unsurprisingly, discrepancies begin to appear for observer angles exceeding $135°$⁠. At such values, the experimental setup gives rise to measurable acoustic reflections and scattering effects behind the propeller.⁴⁴ Nonetheless, excellent agreement between theory and experimental data corresponding to right-handed propellers persists at observer angles above and in front of the HLP. Here too, experimental deviations between the right-handed and left-handed propellers are reported in the wind tunnel measurements. Interestingly, such variations are neither predicted by the surface-vorticity approach outlined above nor in the ANOPP2 computations. A more detailed discussion of this behavior is furnished by Litherland et al.⁴⁴ Apart from this issue, however, Fig. 8 remains convincingly supportive of the validity of the overarching acoustic framework.

To further establish the robustness and sensitivity of the present aeroacoustic formulation in the surface-vorticity solver, three different numerical case studies are undertaken. The first of these focuses on a single five-bladed propeller that is mounted on a contemporary UAM vehicle, as shown in Fig. 9. The UAM vehicle chosen for this example corresponds to the publicly available NASA OpenVSP model of the Joby S4 UAM eVTOL vehicle. The model is transferred from OpenVSP to the solver as a Plot3D geometry with both solid- and thin-surface representations of the propeller blades. The fuselage, wings, and nacelles are Boolean united in OpenVSP before transfer. The propeller is operated at 955 RPM and $α=0°$⁠, while the tip Mach number is set at $Ma=0.4$⁠. Additional propeller parameters are furnished in Table II.

Throughout these simulations, the observer time signal window is taken to be approximately 0.15 s. Moreover, in what concerns this case study as well as the ones to follow, the full two-way aerodynamic coupling is maintained between the overall vehicle and the propeller, as part of the unsteady aerodynamic flow solutions. However, the absence of an acoustic shielding model implies that the acoustic signal emanating from the propeller blades remains unaffected by the shielding effects of the vehicle relative to the observer. In practice, this leads to unrestricted pressure wave propagation throughout the domain. As for the acoustic shielding and scattering effects, they are hoped to be explored in future work.

The acoustic signal histories are generated at 36 microphones that are placed in the form of a circular array at a range radius of 50 m within the plane of the propeller and centered about the propeller hub. Each microphone is separated from its neighbor by an azimuthal angle of $10°$⁠. For each observer, a pressure–time signal is recorded and analyzed to the extent of retrieving the overall sound pressure level (OASPL) at each observer location. Subsequently, the pressure signals obtained at various angular positions are collected and displayed in Fig. 10.

According to this angular distribution, the flow is directed along the $180°→0°$ line into the propeller, which faces $180°$⁠. At first glance, several key characteristics may be inferred from these results. First, the acoustic pressure signal distribution changes azimuthally around the propeller. Second, the minimum signal magnitude is obtained in the front and rear of the propeller, whereas the maximum values appear in the tip plane of the propeller disk. The deviations between these two extremes can be seen in Fig. 10. Note that the numerical error in the acoustic signal is dominant at both $0°$ and $180°$ azimuthal locations, with no physical signal profile visible. At $90°$⁠, strong sinusoidal signals that correlate with the operating RPM of the propeller are observed. These results are consistent with the response of a single propeller in acoustic isolation.¹⁸ 

As described in Sec. II B 1, and further illustrated in Fig. 3, one can rely on so-called acoustic volume sections, which consist of planar arrays of stationary microphones that are distributed around the vehicle, using well-delineated shapes, sizes, and resolutions within the solver. The acoustic field around the vehicle can then be deduced from the resulting contours of sound pressure–time histories obtained across these sections.

For the Joby S4 UAM eVTOL vehicle with a single propeller, an outer spherical radius of 15 m is used to delimit the acoustic field, as shown in Fig. 11. This enables us to visualize the sound pressure distribution around the vehicle while the propeller is operating in the forward-flight mode. For example, the top view section given by Fig. 11(a), which corresponds to an XY plane taken at a distance of 5 m below the vehicle centroid, the sound concentration around the propeller in operation can be clearly seen. However, as the distance from the propeller is doubled to 10 m, the signal is weakened to the extent of becoming barely visible, being masked by inevitable background noise [Fig. 11(b)]. The front view of the vehicle, which is given by Fig. 11(c), clearly shows the circular sound bubble generated by the propeller in the YZ plane taken 3 m upstream of the vehicle, with the strongest signals around the tips of the rotor disk and dissipating radially outwardly. These results are complemented by the XZ section field in Fig. 11(d), where a side view of the aircraft taken 8 m leftward helps to identify the dissimilarities in the signal strength distribution around the vehicle. Given the present cruising speed, the sound generated by the propeller may be seen to trail downstream behind the rotor disk and then dissipate spherically outwardly.

The next case study enables the remaining five propellers on the UAM vehicle to operate, thus illustrating the typical multi-rotor forward-flight DEP and eVTOL industry applications. To do so, all six propellers are operated synchronously in the forward-flight mode at 955 RPM and a tip Mach number of 0.4. The corresponding UAM vehicle configuration is shown in Fig. 12. The surface-vorticity flow solution, which is provided in Fig. 12(b), depicts very distinct vortex stream formations from the interacting propellers. Note that all of the propeller wakes and geometric details are aerodynamically coupled with the fuselage, nacelles, and wings. As such, the wake interference from the front rotor onto the rear rotor can be viewed in Fig. 12(b).

As done in Case 1, 36 stationary microphones are placed in the form of a circular array at three specific radii of 50, 100, and 150 m, respectively; these are placed in a plane 4 m below the base of the vehicle fuselage and centrally with respect to the computed centroid of the vehicle. Each microphone is separated from its neighbor by an azimuthal angle of $10°$⁠. For each observer, a pressure–time signal is recorded and analyzed to retrieve OASPL values at each observer location. These are collected and then displayed in Fig. 13, where several pressure signals measured at 30° intervals are shown. As before, the flow runs parallel to the $180°→0°$ line into the propeller, which faces toward the $180°$ position.

As expected, the loudest signal is observed at right angles to the propeller blades while looking at the tips of the blades. However, the flow signals are now substantially different from those obtained in Case 1 for an acoustically isolated propeller. The net signal at these locations reflects the acoustic in-phase and out-of-phase coupling of the six independent propeller signals. It is noteworthy that the signal in front of the vehicle at $180°$ has a mean-zero signal magnitude with the exception of smaller signal spikes at 0.06 and 0.09 s; these correspond to the in-phase periodic summation of the propeller signals. A more washed-out signal in the opposite direction can be seen at the rear of the vehicle around $0°$⁠. We also note the spikes in the peak signal time positions at the $90°$ location due to the signal travel time to the observer from the six spatially separated propellers. It should be recognized that most of the tonal power stems from the propeller's RPM, with secondary contributions from harmonics related to the blade passage frequency. In this case, as with the first case, the spikes are detected along multiples of the blade passage frequency.

Having examined the acoustic pressure signal at discrete points and the vortex stream formations around the vehicle, attention is now turned to the acoustic field distribution surrounding the vehicle. Similar to Case 1, the actual field can be visualized using acoustic volume sections, as shown in Fig. 14. Of particular note here is that in the XY plane [top view in Fig. 14(a)], the strong signal characteristics are visible in the section located 5 m below the vehicle centroid. By the time the signal reaches 10 m below the vehicle, the signal dissipation becomes quite visible [Fig. 14(b)]. Also visible are the strong signals engendered near the tips of the propeller rotor disks, as seen in the YZ and XZ section fields of Figs. 14(c) and 14(d). Note the phase change of the signal strength upstream and downstream of the rotor disk. Furthermore, and as seen before, an acoustic noise bubble develops around the vehicle near field in the YZ plane, 3 m upstream of the vehicle centroid, where it just touches the nose of the vehicle [Fig. 14(c)]. These visual characteristics are complemented by the side view of the aircraft, which is taken 8 m to the left of the vehicle [Fig. 14(d)].

In addition to pressure field visualization using acoustic sections, the solver has been augmented with the ability to produce cutaways of acoustic spheres. These are showcased in Fig. 15, where two acoustic spheres are generated around the Joby vehicle in the forward-flight mode. In this example, all six 5-bladed propellers are set in synchronous operation, and the images depict two instants of time, particularly, t = 0.04 s and 0.1 s. Moreover, the acoustic pressure values (in Pa) are shown on the surface of the spheres in these visual representations. It should be noted that signal data can be readily computed at the centroid of each spherical surface panel as a function of time. This allows the representation of the complex vehicle signal characteristics on a simplified spherical surface, which can then be post-processed rather straightforwardly for a far-field data analysis. Clearly, through these types of graphical representations, one can achieve a very dynamic visualization of the change in acoustic pressure as a function of time on the surface of the sphere, namely, as the tips of the propeller blades make their closest passes to different areas in the field of study at various observer times.

The last illustrative case is performed on a pusher-propeller DEP UAM vehicle modeled using NASA OpenVSP around the Kittyhawk KH-H1 DEP vehicle depicted in Fig. 16. The purpose here is to demonstrate the differences in acoustic signatures from two highly dissimilar conceptual UAM/DEP designs relative to Cases 1 and 2. In a sense, demonstrating that the surface-vorticity solver can be effectively used to capture these vehicle-level differences may be viewed as a critical requirement for justifying its adoption in a conceptual design.

The KH-H1 vehicle has eight pusher-type propellers in the forward-flight mode shown in Fig. 16. On the one hand, the two pairs of inboard wing propellers are operated at 1430 RPM with a tip Mach number of 0.2. On the other hand, the wingtip propellers are operated at a slightly higher angular speed of 1670 RPM with a tip Mach number of 0.25. As for the canard propellers, they are operated at the highest RPM of 1910, thus leading to a tip Mach number of 0.28. Additional design parameters of the KH-H1 propellers are provided in Table III. In the forward-flight mode, all propellers are allowed to operate synchronously. Therein, the interactions of the slipstream from the front canard propellers with the main wing and nacelles of the wing propellers are quite visible, along with the interference effects of the slipstream with the horizontal stabilizers.

As usual, 36 stationary microphones are placed in the form of a circular array at different range radii of 50, 100, and 150 m, respectively; these are positioned in a virtual plane that is dropped 4 m below the base of the vehicle fuselage and distributed equidistantly from the computed centroid of the vehicle. In this plane, the pressure signals are measured at discrete angular locations that vary from $0°$ to $360°$ in increments of $10°$⁠. These are computed and shown in Fig. 17 at $30°$ intervals.

As before, the flow moves horizontally in a direction from $180°→0°$ into the propeller facing $180°$⁠. The variation in signal intensity and phase identification relative to Cases 1 and 2 is clear. Overall, the signal characteristics that resemble those of Case 2 are realized. However, the presence of eight main propellers in operation compared with six, which actually rotate at three dissimilar blade revolution rates, can be seen to produce a highly jagged acoustic signal at the observer locations. This behavior may be inferred from the wide variety of tones detected. Moreover, the presence of distinct tonal spikes is no longer observed, and this may be attributed to the outcome of signal superposition being prescribed by the phase integrations of various contributing signals. At this juncture, should one be interested in identifying the dense harmonics comprising these signals, a frequency domain analysis will be required, as shown systematically in the companion study by Little et al.⁴ 

For further illustration, acoustic volume sections of the Kittyhawk KH-H1 DEP vehicle with all eight 3-bladed propellers operating synchronously in the forward-flight mode are provided in Fig. 18. Therein, the outer radius of the acoustic volume section is set at 15 m, and the grid resolution is taken to be 40 × 40 in the radial and azimuthal directions, respectively. As before, using a blade pitch angle of $α=5°$ everywhere, the two pairs of inboard wing propellers are operated at 1430 RPM with a tip $Ma=0.2$⁠, the wingtip propellers are operated at 1670 RPM with $Ma=0.25$⁠, and the canard propellers are operated at 1910 RPM with $Ma=0.28$⁠. Figure 18 helps to visualize the acoustic field around the vehicle using the same orthogonal XYZ coordinate system convention. As before, the horizontal X and Y axes run parallel to the fuselage and wings, respectively, and Z represents the vertical direction. In the top view given by Fig. 18(a), the signal characteristics are captured in a plane situated at 5 m below the vehicle. In Fig. 18(b), one may note a strong signal attenuation while moving to a plane situated 10 m below the vehicle. As usual, an acoustic noise bubble may be observed around the vehicle in the YZ plane, 3 m upstream of the vehicle centroid [Fig. 18(c)]. These visual aids are corroborated by the side view of the aircraft, which is taken at 8 m left of the aircraft, as shown in Fig. 18(d). They are also consistent with the three-dimensional acoustic spheres produced around this vehicle and provided in Fig. 19 at two different time steps, namely, at t = 0.04 and 0.2 s, respectively.

It may be instructive to note that the acoustic spheres around the KH-H1 vehicle have a relatively short radius, thus allowing the acoustic signal to quickly reach steady-state conditions on the spherical boundary, as shown in Fig. 19. The presence of aft-facing DEP propellers on the main wing seems to have a strong bearing on the acoustic pressure contours. These exhibit higher pressure magnitudes aft of the aircraft and negative pressures upstream of the propellers. Clearly, the presence of multiple propellers tends to promote a uniform acoustic pressure field wherein the blade-induced perturbations of each propeller are effectively suppressed. The resulting behavior is also suggestive that various DEP configurations can lead to dissimilar acoustic additive and dissipative interference characteristics that can be, at least in principle, effective at achieving acoustic signal reduction.

This study details the effort to integrate and illustrate the application of the aeroacoustic Farassat Formulation 1A within a robust surface-vorticity solver. The combination of a fast and user-friendly flow solver with an aeroacoustic toolbox has led to a capability in which UAM and UAV vehicles can be readily imported from a large variety of contemporary engineering tools. Given the high efficiency of the overarching framework, accurate solutions obtained in a matter of minutes are shown to include not only the flow characteristics, vorticity streaks, and loads all around the aircraft but also the acoustic signatures at essentially any observer point of interest. As part of this effort, extensive post-processing tools have been systematically conceived and developed to aid the designer in meeting acoustic signature requirements. The underlying capability is further validated through comparisons with recent numerical simulations and experimental measurements of NASA's X-57 propeller, thus showcasing reasonable agreement across a wide range of observer angles and acoustic pressure levels. In all cases considered so far, the present simulation results are found to agree reasonably well with openly available experimental measurements and computations obtained at subsonic speeds. In this process, however, some known limitations of the F1A formulation are identified; these include its inability to resolve the high-frequency broadband noise or capture the nonlinear quadrupole noise terms that become appreciable at high Mach numbers. In addition to a preliminary study of NASA's RVLT Tiltwing UAM eVTOL vehicle with eight 5-bladed propellers, three specific UAM representative vehicles are selected to showcase the multi-rotor acoustic signal prediction capabilities. This is accomplished using extensive acoustic pressure characterization studies of the Joby S4 UAM eVTOL concept with either one or six synchronously operating five-bladed propellers as well as the Kittyhawk KH-H1 DEP vehicle with eight 3-bladed propellers.

Since this effort focuses on the strategic and programmatic implementation of the F1A formulation, the frequency-domain decomposition, post-processing, and evaluation of several certification-driven acoustic metrics of interest to the Federal Aviation Administration (FAA) and International Civil Aviation Administration (ICAO) are provided separately in the companion article by Little et al.⁴ As for the effects of transition and hover flight phases of eVTOL-class aircraft, they are hoped to be addressed in future analysis. Work in this direction is actually underway, namely, through the NASA Advanced Air Mobility program. As part of this ongoing effort, the acoustic field-test results for transition, overflight, and hover conditions will be extensively used to further refine the framework under development.

The authors wish to acknowledge the cooperation and support of the NASA researchers and the NASA SBIR Office on this effort. The primary work for this study was performed within the context of the NASA STTR 80NSSC20C0586 Phase I contract activity.

The authors have no conflicts to disclose.

Vivek Ahuja: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Daniel S. Little: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Resources (supporting); Software (equal); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Joseph Majdalani: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (supporting); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (equal). Roy Hartfield: Conceptualization (supporting); Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Methodology (supporting); Project administration (equal); Resources (equal); Software (supporting); Supervision (equal); Validation (supporting); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.