Acoustic analogy methods are used as post-processing tools to predict aerodynamically generated sound from numerical solutions of unsteady flow. The Ffowcs Williams–Hawkings (FW–H) equation and related formulations, such as Farassat’s Formulations 1 and 1A, are among the commonly used analogies because of their relative low computation cost and their robustness. These formulations assume the propagation of sound waves in a medium at rest. The present paper describes a surface integral formulation based on the convective wave equation, which takes into account the presence of a mean flow. The formulation was derived to be easy to implement as a numerical post-processing tool for computational fluid dynamics codes. The new formulation constitutes one possible extension of Farassat’s Formulation 1 and 1A based on the convective form of the FW–H equation.
As a branch of computational fluid dynamics (CFD), computational aeroacoustics (CAA) is dedicated to the study of aerodynamically generated sound. One common approach in CAA is to perform numerical simulations of the flow within a limited region of interest where sources are located, i.e. the nearfield. The sound radiated to the farfield is then calculated by solving a wave equation. This indirect approach, which is referred to as the acoustic analogy method in the literature, is often based on the assumption of linear propagation of waves from the nearfield to the farfield.
Depending on the motion of the observer and the ambient fluid speed, three possible problem types may be encountered, as summarized in table 1. The first type, referred to as fly-over, involves sound perceived by a stationary observer in a stationary ambient medium. The source may be either stationary or moving (e.g. the flight of an aircraft over a microphone array). The second type, referred to as moving-observer, deals with problems in which the observer is moving (e.g. microphone mounted on the fuselage of a helicopter or an aircraft). In the third type, referred to as wind tunnel, sound is received by a stationary observer in a moving medium such as sound measurements in wind tunnels.
A Galilean transformation may allow a problem of one type to be transformed into an equivalent problem of a different type. For example, a wind-tunnel problem in which the medium is moving at velocity U0 may be transformed into a moving-observer problem, where the observer and the source are translating at velocity −U0 with the medium being stationary. Such transformations are aimed at simplifying the problem and reducing the computation cost.
The Ffowcs Williams–Hawkings (FW–H) equation and related formulations, such as Farassat’s Formulation 1 and 1A (Farassat 1975; Farassat & Succi 1980; Brentner 1997a), feature a relatively low computation cost, and they are fairly robust. These formulations have been successfully used for a wide range of fly-over and moving-observer problems (cf. Brentner & Farassat 2003; Lyrintzis 2003).
Formulations 1 and 1A do not explicitly take into account the presence of a mean flow for wind-tunnel problems. One common practice is to transform a given wind-tunnel problem into a moving-observer problem where the observer is assumed to be moving at a constant speed in a quiescent environment (Brentner & Farassat 2003). One alternative approach is to explicitly take the presence of a mean flow into account by solving a convective wave equation. Besides its mathematical elegance, such formulation helps in distinguishing mean flow effects on from sound generation phenomena.
The objective of the present paper is to present an extension of the classical FW–H equation and its subsequent formulations (Farassat’s Formulation 1 and 1A) to wind-tunnel problems based on the solution of the convective wave equation. As such, it is referred to as Formulation 1C. The formulation is derived to be simple to implement and easy to parallelize. Formulation 1C is specifically designed to make efficient use of the nearfield data obtained from a CFD code.
The rest of the paper is organized as follows. The FW–H formulation for fly-over problems and the derivation of its convective form for wind-tunnel problems are reviewed in §2. The derivation of the new formulation is presented in §3. The numerical implementation and verification cases are discussed in §4.
2. Theoretical background
(a) Acoustic analogy methods for fly-over problems
The concept of the acoustic analogy was pioneered by Lighthill (1952) in his classical study of the sound generated by a turbulent jet (Lighthill 1954). Lighthill’s idea was to recast the exact equations of motion, i.e. the conservation of mass and momentum, as an inhomogeneous wave equation. The problem of calculating the turbulence-generated sound is then equivalent to solving the radiation of a distribution of sources into an ideal fluid at rest.
Lighthill’s equation does not take into account the presence of solid surfaces in the field. A more general formulation which is applicable to the cases where a stationary solid surface is present in the flow field was developed by Curle (1955), and later generalized by Ffowcs Williams & Hawkings (1969) for moving and permeable surfaces. Di Francescantonio (1997) and Brentner & Farassat (1998) demonstrated the advantages of the permeable FW-H surface in CAA. The relationship between the permeable FW–H equation and Kirchoff methods (Farassat & Myers 1988; Lyrintzis 1994) was highlighted in Brentner & Farassat (1998).
The FW–H acoustic analogy involves enclosing the sound sources with a control surface that is mathematically represented by a function, f(x,t)=0. The acoustic signature at any observer position can be obtained from the FW–H equation2.1where [ ]τe denotes evaluation at the emission time τe. The source terms under the integral sign are2.2 2.3and Tij is referred to as Lighthill’s stress tensor,2.4The vectors u and v are the flow and the surface velocities, respectively. See table 2 for the definition of other quantities. The compression tensor P′ij is defined as2.5The third term integral should be evaluated in the region outside of the surface, f>0, to account for the contribution of sound sources outside the FW–H surface. If the FW–H surface coincides with stationary solid boundaries, the flow velocity, u, is equal to the surface velocity v, which reduces the FW–H equation to Curle’s equation (cf. Howe 1998).
The presence of both temporal and spatial derivatives with respect to the observer frame of reference may raise some problems for the numerical implementation of the FW–H equation in the form of equation (2.1). The FW–H equation was later revisited by Farassat (1975), Farassat & Succi (1980) and Brentner (1986, 1997b) and new formulations that are better suited for numerical implementations were introduced. These formulations, known as Formulations 1 and 1A, are widely employed in rotor and propeller noise studies.
An important point to note in equation (2.1) is the distinction between the emission time and the observer time. The left-hand side (LHS) of equation (2.1) denotes the pressure fluctuation perceived by the observer, located at x, at time t. Equation (2.1) implies that the observer pressure fluctuation is obtained through the superposition of contribution from spatially distributed sources denoted by right-hand side (RHS) terms. Elemental contributions from each source are obtained by evaluating the RHS at the emission time, τe, which is different from t. Different sources may also have different emission time values τe depending on their position relative to the observer.
(b) Retarded time versus advanced time algorithms
Two approaches may be followed for the numerical implementation of equation (2.1) based on the distinction between the emission and observer times. In the retarded time approach (Brentner 1986, 1997b; Brentner & Farassat 2003), the solution is evolved in the observer time frame. In other words, the observer time t and location x are kept constant at each time step. The emission time is calculated for each source location, y. The associated flow properties are obtained using a temporal interpolation to obtain the required information at the emission time τe. This approach requires uploading the flow-field variables of multiple simulation snapshots per each observer time step calculation. Thus, the retarded time approach involves a large memory usage and computational cost associated with the interpolations of the time-resolved CFD results.
In the advanced time approach (Leishman 1996; Lyrintzis & Xue 1997; Casalino 2003; Lyrintzis 2003), also known as the source time-dominant algorithm (Brentner & Farassat 2003), the source time is chosen as the primary variable to evolve the solution. At each flow-field time step, matrices of instantaneous flow properties from the CFD solution are loaded into memory. The RHS of equation (2.1) is then evaluated at each point in space along with the arrival time value (advanced time). The resulting contribution is then added to the time history of the pressure received by the observer at the proper advanced time. Once the contributions from all points in the flow field have been calculated, the flow properties at the next source time step are loaded. As the calculation is evolved along the source time axis, the observer’s sound pressure history is constructed. The advanced time approach is better suited than the retarded time, when large datasets provided by CFD codes are used as input data. The source time algorithm results in significant reduction in memory requirement, and is intrinsically parallelizable (Brentner & Farassat 2003; Lyrintzis 2003). A comparison between the cost of the two approaches was reported by Brès et al. (2004).
(c) Wind-tunnel problems
The classical FW–H equation—equation (2.1)—and subsequent formulations 1 and 1A do not explicitly take into account the presence of a mean flow for wind-tunnel problems. One solution is to transform the given wind-tunnel problem into a moving-observer problem, where the observer is assumed to be moving at a constant speed in a quiescent environment. One alternative approach is to explicitly take the presence of a mean flow into account by solving a convective wave equation (cf. Morino 1974, 1985). The second approach explicitly takes the mean flow into account in a cohesive mathematical framework. The resulting formulation provides some physical insight on the mean flow effects on the propagation of acoustic waves, and the presumed sound generation mechanism. It may also help to relate the results from wind-tunnel measurements with those obtained from fly-over experiments. The formulation also allows the clarification of some derivation details. In particular, also working with the advanced time algorithm, Casalino (2003) suggested to change the temporal derivatives with respect to the observer time of the non-convective FW–H equation into a Lagrangian derivative:2.6where x1 is the direction of the mean flow. This seems to be in contradiction with Brentner & Farassat (2003) on using the Formulations 1 and 1A for the moving-observer problems. The formulation of the present paper circumvents this difficulty by directly taking into account the presence of a mean flow rather than solving a moving-observer problem.
(d) The convective FW–H equation
The derivation of the convective FW–H equation was first introduced by Wells & Han (1995). The original derivation of Wells & Han (1995) did not include the quadrupole noise term, and was for a moving-observer problem by considering a moving frame of reference. This situation is equivalent to a stationary observer in a moving medium, i.e. wind-tunnel problems. For the sake of completeness and to clarify the notation, the convective FW–H equation is directly derived for a wind-tunnel problem in this section. Moreover, the quadrupole noise term is retained.
Consider the motion of an acoustic perturbation in a medium moving at constant velocity U0. The flow velocity at each point is U0+u,1 where u is the local perturbation velocity. The continuity equation,2.7can be simplified to2.8where the primed variables denote perturbed quantities. The conservation of momentum,2.9is recast as2.10The closed surface is mathematically represented by f(x,t)=0 (cf. figure 1).2 This surface may be stationary or moving.
Without loss of generality, f may always be defined such that f(x,t)>0 represents the region outside the surface, f(x,t)<0 represents the region enclosed by the surface and |∇f|=1 on the surface. The latter assumption is not necessary, but makes the derivations easier (Farassat 1994; Brentner & Farassat 2003).
The continuity equation in the region outside the surface is then2.11where H is the Heaviside function. Moving H(f) inside the differential operators yields32.12Using2.13and2.14equation (2.12) is recast as2.15where2.16Similarly, the conservation of momentum is written as2.17where2.18is Lighthill’s stress tensor and2.19In equation (2.18), the stress tensor Pij is defined as2.20
Differentiation of equation (2.15) with respect to time, calculation of the divergence of equation (2.17) and subtraction of the latter from the former yields2.21Finally, equation (2.15) is used to replace the second term on the LHS of equation (2.21) and obtain the convective FW–H equation,2.22The convective FW–H equation is an inhomogeneous convective wave equation similar to the well-known, classical, non-convective FW–H equation. The first two terms on the RHS of equation (2.22) are a (convective) monopole term, also known as thickness source, and a dipole term which is also called the loading source. The last term is a quadruple source term, which is typically small when compared with the other contributions. This term is also more challenging to compute and is often neglected. The specific form of the convective wave operator2.23takes into account the presence of a mean flow and reduces to the simple wave operator when U0=0 as in the classical FW–H equation. The solution of equation (2.22) requires a Green’s function that takes the presence of a mean flow into account, i.e. a convective Green’s function. The thickness, Qj, and loading, Lij tensors include terms with the mean flow velocity and are slightly different from their counterparts in the non-convective FW–H equation. Lighthill’s tensor definition retains the same as for the classical non-convective formulation.
Without loss of generality, it can be assumed that the mean flow velocity is along the positive x1-direction.4 Assuming a subsonic mean flow, the three-dimensional free-space Green’s function for the convective wave equation is (cf. eqns (1.91) and (1.92) of Blokhintsev 1956)2.24where2.25 2.26and2.27Using the above Green’s function yields the solution to the convective FW–H differential equation,2.28where2.29Equation (2.28) shows that the thickness of noise temporal derivative in the non-convective FW–H equation is indeed replaced by a Lagrangian derivative as suggested by Casalino (2003); however, this substitution is not sufficient to obtain correct results. In addition, Qj and Lij should be modified to include U0j terms, as shown in equations (2.16) and (2.19). Moreover, the convective Green’s function should be used, and not the free-space Green’s function.
3. Formulation 1C
As for the classical non-convective FW–H equation, the numerical solution of the convective FW–H equation (2.28) is challenging because of the concurrence of spatial and temporal derivatives. This problem may be circumvented as presented next. The contribution of quadrupole noise terms is neglected in the following,5 and the thickness and loading terms are treated separately.
(a) Thickness noise
The thickness noise contribution, p′T, is obtained from3.1The surface f=0 may be defined in a frame of reference fixed to the surface, η. Any point on the surface is defined by a fixed value of the variable, η, irrespective of the motion of the surface. The motion of each point on the surface is then fully described in terms of the translation and the rotation of the η frame of reference. It is also assumed that the FW–H surface does not undergo deformation, expansion or contraction. Since the transformations3.2and3.3are isometric, the Jacobian of the transformation is unity. Hence, equation (3.1) may be written as:3.4To obtain a formulation that is suitable for numerical implementation, the spatial derivative ∂/∂x1 must be converted into a temporal derivative. All terms in the integral are functions of η and τ only, with the exception of δ(g)/R*, which depends on x and t as well. The combination of3.5and3.6yields3.7where the radiation vector is3.8and3.9These relations may be used to recast equation (3.1) as3.10The next step is to change the variable from τ to g. By definition, ∂g/∂τ is obtained from3.11where MR is defined as3.12Using equation (3.11), we have3.13Hence, equation (3.10) is simplified to3.14where […]ret denotes the evaluation at the retarded (emission) time3.15
The quantity R is no longer the geometric distance between the observer and the source, but the acoustic distance between the two, as defined in equation (2.25) (Wells & Han 1995). In the limit of zero mean flow velocity, the acoustic distance and the geometric distance are equal. The observer temporal derivative, ∂/∂t, may be evaluated numerically using a backward difference in time.
Some difficulties may arise in the reconstruction of the signal using the advanced time approach. The input flow properties are usually provided at equally spaced time steps in the source time domain; however, the uniform discretization of the source time domain does not necessarily yield a uniform discretization of the observer time domain because of the Doppler effect. This inconvenience may be avoided by moving the observer temporal derivative inside the integral, and transforming it into a source temporal derivative (cf. appendix A for details of the derivation).
The thickness noise contribution to the farfield sound pressure is then equal to3.16where dots over quantities denote temporal derivatives with respect to the source time. The term is obtained numerically, while other temporal derivatives are obtained either numerically or analytically.
(b) Loading noise
The loading noise is obtained from3.17As for the thickness noise, all terms under the integral are functions of η and τ only, with the exception of δ(g)/R*. Therefore, equation (3.7) is used to transform the spatial derivative, ∂/∂xi, into the observer temporal derivative, which yields3.18Equation (3.13) is used to further simplify the above equation to3.19It can be verified that letting U0=0 recovers the loading noise term in Farassat’s Formulation 1. The loading noise contribution may also be written in terms of source temporal derivatives as follows:3.20This equation is equivalent to the loading noise term of Farassat’s formulation 1A, and is convenient for the implementation of the advance time algorithm.
(c) The special case of ‘wind tunnel’
In the particular case where both the source and the observer are stationary in a wind tunnel, simplifications in Formulation 1C lead to increased computational efficiency.
In Formulation 1C, the distances R and R*, defined by equations (2.25) and (2.26), are constant and do not vary with time. The same observation applies to the radiation vector components , the local normal vector components ni, and MR=0. Therefore, these variables can be evaluated and stored at preprocessing, rather than being computed at every time step, and their source temporal derivative is zero. Likewise, the source Mach number M and the unit normals n are not functions of time, so and . These simplifications lead to the following expressions for the thickness noise3.21and the loading noise3.22The terms containing R*2 vanish quickly in the farfield as , and thus they are significant only in the nearfield. It is important to note that the amplitude of the thickness term contribution to the farfield noise decreases with increasing mean flow Mach number. Therefore, loading noise becomes the dominant source term for high Mach number flows.
4. Numerical implementation and verification
Formulation 1C was implemented in an object-oriented code written in C++. The input to the code is a surface mesh, which defines the FW–H surface, and the time-accurate flow solution obtained from a distinct CFD code. Surface elements (panels) may be triangular, quadrilateral or polygonal. The flow properties are specified at the centre of each panel (for first-order accuracy) or at the vertices of the panel (for higher order accuracy) at each time step of the solution. The temporal derivatives are obtained from a finite difference approximation. Finite-element type shape functions are used for the spatial interpolation. The integral terms are then obtained from a Gaussian quadrature with an appropriate number of points to obtain the desired order of accuracy. The FW–H code was written based on the advanced time algorithm. The canonical problems of sound radiation by a monopole and a dipole in a uniform mean flow were solved for the numerical verification of the formulation. Source terms were computed over a fictitious, closed surface surrounding the source.
In all test cases, the mesh resolution of the FW–H surface and the source and observer time steps are chosen such that the wavelength and period of the radiated waves are resolved with at least eight points.
(a) Test case 1. Stationary monopole in a moving medium
The monopole sound-field is characterized by a simple harmonic velocity potential function (Lockard 2002)4.1The induced particle velocity is obtained from4.2The induced pressure and density, for the case with a mean flow U0 along the x1-direction, are related to the velocity potential through4.3and4.4respectively.
To keep a focus on the formulation and its implementation, and to avoid any bias related to the accuracy of the flow solver, the flow properties on the FW–H surface were obtained from the exact solution of the flow field generated by the monopole, i.e. equations (4.1)–(4.4). These properties were then used as the input to the code.
The monopole was located at the origin of the coordinate system. The velocity potential amplitude was A=1 m2 s−1. The pulsation frequency was 5 Hz. The radiated sound pressure was calculated at a geometrical distance of 340l from the source, where the characteristic length of the FW–H surface, l, was set to unity. The speed of sound was 340 m s−1. Figure 2 shows the directivity pattern of the sound pressure measured in the presence of a uniform mean flow along the the positive x1-direction for M0=0.5 and M0=0.85. In this polar plot, the radius indicates the RMS pressure at the observer location, while the polar angle corresponds to the geometrical angle (measured from the x1 axis) between the source and the observer. In a moving medium, the directivity pattern is directional owing to the convective amplification caused by the mean flow. Convective amplification causes the sound pressure level upstream of the source to be greater than that downstream at the same distance from the source. This convective effect increases with the mean flow Mach number. The ratio of the sound pressure at θ=180° (upstream) and that at θ=0° (downstream) is equal to the theoretical value of (1+M0)/(1−M0) with less than 0.3 per cent error.
The FW–H equation applies in principle for observer locations in both the nearfield or the farfield, and involves no farfield approximation. The accuracy of the sound prediction in the nearfield was also verified. Figure 3 shows the directivity pattern measured in the nearfield with microphones located at r=5l. The exact solution and the numerical predictions are in very good agreement. The small deviation from the exact solution is due to the fact that measurements were made very close to the FW–H surface, where the size of the panels (mesh cells) is comparable to the distance between the source and the observer. This means that the panels are no longer acoustically compact. A finer surface mesh yielded more accurate results (results not shown here).
(b) Test case 2. Stationary dipole in a moving medium
The second validation test case is sound radiation from a point dipole in a moving medium. The dipole axis was aligned with the x2-axis. The velocity potential for such dipole in a mean flow is4.5The particle velocities and pressure fluctuations were obtained from equations (4.2) and (4.3), respectively. The same potential amplitude, frequency and position as test case 1 were prescribed. Figure 4 shows the directivity pattern of the sound pressure at r=30l. The mean flow causes the direction corresponding to the maximum acoustic pressure to move in the upstream direction. The FW–H code and the exact solution yield nearly identical results.
(c) Test case 3. Rotating monopole in a moving medium
The last test case is a rotating monopole radiating in a moving medium as shown schematically in figure 5. This case was designed to validate the accuracy of Formulation 1C to predict the general case of sound radiated by moving sources in uniformly moving media. This case features radiation effects similar to those of the thickness source term of a fan or a helicopter rotor blade measured in wind tunnels. A monopole was placed initially at (0.7l,0,0), surrounded by an FW–H surface moving with the source. The monopole rotated around the x3-axis with an angular speed of 2π (rad s−1). The potential amplitude was A=1 m2 s−1 and the pulsation frequency was 5 Hz as for the previous test cases. The ambient medium had a Mach number of M0=0.5 moving in the positive x1-direction. Figure 6a shows the time history of the sound pressure perceived by an observer located at (2l,0,0). Excellent agreement between the exact solution and the farfield prediction confirms the validity of Formulation 1C and its implementation.
The effect of the mean flow is illustrated in figure 6b, where the time histories of the sound pressure measured upstream (θ=0°) and downstream (θ=180°) of the monopole at a distance of r=2 l are compared. As expected, the amplitude of the peak pressure for the upstream observer is greater than that for the downstream observer because of convective amplification. The pressure ratio is no longer (1+M0)/(1−M0), as in case 2, because of the motion of the monopole.
An acoustic analogy formulation based on the convective form of the FW-H equation was introduced. The formulation, called Formulation 1C, was used to simulate wind-tunnel acoustic measurements of stationary or moving sources. The same formulation is also useful for moving-observer problems.
In comparison with the classical non-convective FW–H equation, the temporal derivative of the thickness noise changes to a Lagrangian temporal derivative in the convective FW–H formulation. The monopole (thickness) and dipole (loading) sources also change slightly. The quadrupole term retains the same form as in the classical FW–H equation. The convective FW–H equation was solved using the Green’s function suggested by Blokhintsev (1956).
Formulation 1C is amenable to the advanced time approach implementation, which may significantly reduce storage and computation requirements. The validity of formulation 1C and its implementation was demonstrated for the canonical test cases of a monopole and a dipole radiating in a moving medium. A test case for a rotating monopole was also investigated. In all cases, numerical results were in very good agreement with the exact solutions.
This work was sponsored by Exa Corporation. The authors would like to thank Drs David Freed, Frank Perot, Richard Shock from Exa Corporation and Phoi-Tack Lew from McGill University for their help and comments.
Appendix A. Derivation of Formulation 1C in terms of source temporal derivatives
The details of the derivation of the thickness noise in terms of the source temporal derivative, namely equation (3.16), are shown. The derivation of the loading noise term, equation (3.20), follows the same steps, and is not presented.
Consider a general expression of the formA1Taking the temporal derivative inside the integral, and using the Leibniz rule, yieldsA2The second integral on the RHS of the above equation is zero. This can be proven as follows. Consider the spherical coordinate system (R*,ϕ,θ), whereA3andA4The above transformation can be used to replace y and R in terms of R*A5The second integral can then be recast asA6The above equation can be further simplified toA7Since R*δ(R*−ϵ)=ϵδ(ϵ−R*), the above integral vanishes as . HenceA8The above identity means that we can take the temporal derivative inside the integral. The integrals of equation (3.10) are then contracted to obtain equation (3.14), and the differential operator is moved inside the integral to obtainA9The last step is to transform the differentiation with respect to observer time into a differentiation with respect to source time,A10where subscript x means that x is kept constant in the above differentiation. Since emission time τe is the root of g(x,t;y,τ)=0,A11which yieldsA12Using the above identity, the final form of the thickness noise can be deduced as presented in equation (3.16). Similarly, the loading noise takes the form that is presented in equation (3.20).
↵1 The bold font is used to denote vector quantities. Such quantities are also denoted as tensors of rank one and follow Einstein’s index notation.
↵3 Differentiation across a discontinuous function involves the use of generalized functions, also known as distributions in functional analysis and partial differential equations. See, for example, Renardy & Rogers (2004) or Farassat (1994, 2000) for an introduction to generalized functions.
↵4 If the mean flow is not along the x1 direction, the reference frame can be rotated to satisfy this condition.
↵5 Di Francescantonio (1997) and Morgans et al. (2005) discuss in detail that the permeable control surface can be chosen, such that the important noise sources are enclosed without the need to calculate the quadrupole noise.
- Received March 30, 2010.
- Accepted May 26, 2010.
- © 2010 The Royal Society