## Abstract

Noise generation in a subsonic round jet is studied by a simplified model, in which nonlinear interactions of spatially evolving instability modes lead to the radiation of sound. The spatial mode evolution is computed using linear parabolized stability equations. Nonlinear interactions are found on a mode-by-mode basis and the sound radiation characteristics are determined by solution of the Lilley–Goldstein equation. Since mode interactions are computed explicitly, it is possible to find their relative importance for sound radiation. The method is applied to a single stream jet for which experimental data are available. The model gives Strouhal numbers of 0.45 for the most amplified waves in the jet and 0.19 for the dominant sound radiation. While in near field axisymmetric and the first azimuthal modes are both important, far-field sound is predominantly axisymmetric. These results are in close correspondence with experiment, suggesting that the simplified model is capturing at least some of the important mechanisms of subsonic jet noise.

## 1. Introduction

Research in the field of jet noise has recently been dominated by the development and application of numerical approaches including direct and large-eddy simulation. Examples include the direct numerical simulation (DNS) of Freund (2001), and the large-eddy simulations of Bogey *et al*. (2003) and Shur *et al*. (2005). This work has demonstrated that jet noise is increasingly tractable as a solution of the compressible Navier–Stokes equations. For reasons of computational expense, the simulations are limited with respect to resolution of the nozzle lip (which may be important for flow control) and in the frequency range of the simulated spectra. However, significant progress has been made, with the best simulations now able to predict noise to within 2–3 dB (Shur *et al*. 2005) for jet Strouhal numbers up to 1.5. Superficially, one might suppose that there is now little place for theoretical work on this problem and that, in the context of the present journal issue on the nature of fluid dynamics research in the twenty-first century, many other problems may similarly be considered ‘solved’ when the Navier–Stokes equations are simulated numerically. In this paper, we aim to show that, even for a problem in which DNS has been successfully carried out, there is still an important role for theory in providing physical insight (which is rather more satisfying than contemplating terabytes, or in the not too distant future petabytes, of data), in interpreting the numerical results and constructing less expensive predictive models. Indeed, the present line of research (cf. Sandham *et al*. 2006*a*) started when we tried to interpret results from a DNS of a subsonic jet. We then ran simplified simulations and theoretical models to try to understand the flow physics. Along with traditional methods of enquiry, we speculate that this type of study will become more common as numerical simulation codes and databases become more widely available in the applied mathematics community.

In contrast to supersonic jets, where the primary flow instability directly radiates sound (Tam & Burton 1984; Wu 2005), subsonic jet noise is less amenable to theory (Tam & Morris 1980). However, there is a body of experimental work that suggests that instability modes (solutions of the compressible Orr–Sommerfeld equation for example) are present in jets, and thus the possibility exists that they can be exploited in prediction schemes, as discussed in Huerre & Crighton (1983). Laufer & Yen (1983) located the source of sound with the nonlinear saturation of instability waves and noted that, while the near-field pressure varied linearly with disturbance amplitude, the acoustic far field had a quadratic dependence. Stromberg *et al*. (1980) made experiments in jet Mach number *M*=0.9 and *Re*=3600 based on jet diameter and also suggested that their data supported a nonlinear mechanism of noise generation. This experiment is attractive for analysis and simulation, not least because the authors included velocity profiles from the region near the nozzle exit. The same configuration was the subject of a DNS study by Freund (2001) in which good agreement was found with the experiment. The simulation added more complete flow field data and we choose this case for the present study. Although the role of instability waves is naturally clearest in work at lower Reynolds numbers, instability waves have been detected in turbulent jet experiments by Suzuki & Colonius (2006) at Reynolds numbers up to 10^{6}, suggesting that the basic mechanisms may be inviscid in nature.

A widely used model for sound generation in shear flows is the Lilley–Goldstein analogy (Lilley 1974; Goldstein 2001), which extends Lighthill's approach (Lighthill 1952; Morfey & Wright 2007) to be applicable to disturbances on a parallel shear flow. The Lilley–Goldstein equation can be written in cylindrical polar coordinates with radial direction *r* and axial direction *z* as follows:(1.1)where , *p* is the wave variable (denoted by *π* in Goldstein 2001) and the parallel shear flow is defined by the axial velocity , density and sound speed . At the highest order for nearly isothermal jets, the source terms are given by(1.2)Components of the source term in cylindrical polar coordinates can be found, for example, in Tester & Morfey (1976).

In Sandham *et al*. (2006*a*) a model problem was used to show how products of linear instability modes could be used as source terms, leading to potentially useful predictions of the main characteristics of the sound field. Considering first the temporal case, with real wavenumber and complex frequency, any quadratic mode interaction can be split into a term involving the sum of the wavenumbers and another involving the differences. The difference in wavenumber of two instability modes gives a longer wavelength of the radiated sound field compared with the primary instability wavelength, in accordance with the numerical solutions of the compressible Navier–Stokes equations. Figure 1 sketches how this idea is applied to the spatially evolving single-stream round jet considered in this paper. In this case, the frequencies are real and the wavenumbers complex. Two spatially developing unstable waves inside the jet with different frequencies *f*_{1} and *f*_{2} may interact via the source term in (1.1) to give effective sound radiation with frequency *f*_{1}−*f*_{2}. Sandham *et al*. (2006*b*) showed that wave packets based on growing and decaying envelope functions are inefficient at radiating sound for low Mach numbers, slow saturation time-scales and low frequencies; hence interactions involving the summed frequencies will be less effective radiators than those involving difference frequencies.

In this paper, we extend the method from Sandham *et al*. (2006*a*) to spatially developing waves. Wavenumbers and growth factors of instability modes are obtained from the compressible parabolized stability equations (PSE) with viscous and non-parallel terms included. Nonlinear interactions involving different frequencies and azimuthal wavenumbers are used to drive a wave equation to determine the radiation effectiveness of the various possible interactions.

## 2. Formulation

### (a) Instability modes and interactions

In cylindrical coordinates with *z* the streamwise direction, *r* the radial direction and *θ* the azimuthal direction, we consider the interaction of two spatially developing instability modes *u*_{jk} and *v*_{lm} with different (real) frequencies *ω*_{j} and *ω*_{l} and azimuthal wavenumbers *n*_{k} and *n*_{m}. Lengths are scaled with the jet diameter, velocities with the jet core velocity, and density *ρ* and temperature *T* with their respective jet core values. The jet core is defined upstream of the nozzle. The indices *j*, *k*, *l* and *m* will be varied to cover all possible quadratic interactions. The modes evolve according to(2.1)and(2.2)Here c.c. denotes complex conjugate, and are complex eigenfunctions and the complex wavenumber has been split into a real part *α* and an imaginary part −*σ*, both of which may vary in the *z* direction. In later PSE calculations a *z* dependence will be included in the mode shape terms; but the interaction approach can be developed by assuming traditional Orr–Sommerfeld eigensolutions. By direct multiplication, it is possible to write down the source terms for the wave equation (1.1), for example the axial quadrupole term(2.3)where ^{*} denotes complex conjugate and *A*^{+} and *A*^{−} are sum and difference mode amplitudes, given respectively by(2.4)and(2.5)It can be seen that in both cases the growth rates *σ* add, while the wavenumbers *α* can either add or subtract. A basic radiation pattern for any mode interaction can be found by solving a wave equation forced by (2.3). A similar procedure can be carried out for all the other source terms arising from nonlinear interactions, which may involve derivatives of the eigenfunctions.

### (b) Parabolized stability equations

In the present paper, we determine *α*, *σ* and by solving the PSE. This is an efficient method of generating mode shapes and wavenumbers in convectively unstable flows. With PSE it is also straightforward to incorporate non-parallel and disturbance evolution effects that would not be included in the simpler linear stability approach based on parallel flow. The latter approach has also been tested and has yielded similar results, suggesting that the additional effects are not critical to the mechanisms that are described here. The PSE system of equations (Herbert 1997) can be written in concise form as follows:(2.6)where ** K**,

**and**

*L***are matrices and is a vector of mode shapes. The matrix**

*M***contains derivatives in the**

*L**r*direction, which are computed using a sixth-order compact difference scheme on a stretched mesh, extending to

*r*=

*L*

_{r}in the radial direction. The contents of the matrices can be found, for example, in Hein (2005). PSE has recently been applied to jet noise by Cheung

*et al*. (2007): linear PSE was found to be able to predict noise from supersonic jets, while nonlinear PSE, together with an acoustic analogy, was needed for subsonic jet noise. In the present work, by contrast, we use only linear PSE and compute the nonlinear interaction terms separately.

The solution starts at *z*=0 with prescribed , *σ* and *α* from the compressible Orr–Sommerfeld equation for a given *ω*. The Orr–Sommerfeld problem is solved with the same compact finite-difference scheme in *r* and a direct matrix eigenvalue method to find all the eigenvalues and eigenvectors, from which the most unstable mode is selected. After premultiplying by *M*^{−1}, equation (2.6) is advanced in *z* using first-order backward differences. At each step the wavenumber *α*+i*σ* is adjusted iteratively, so that(2.7)This removes the ambiguity of *z* dependence in both the mode shape term and *α* and keeps most of the growth in *σ*.

It is well known that PSE suffers from a limitation on streamwise step size, whereby very small step sizes lead to instability. Various methods have been prescribed to improve this (see Schmid & Henningson 2001). In the present work, PSE was found to be able to cope with smaller streamwise step sizes Δ*z* (without being detrimental to the accuracy) if the streamwise pressure gradient term is omitted from ** M**, as suggested in Herbert (1994). The PSE method was validated against linear stability theory for parallel flows and has also been used in a study of attached boundary-layer flows by Yao

*et al*. (2007).

### (c) Solution of the Lilley–Goldstein equation

The acoustic radiation from the various mode interactions can be found by solving the wave equation (1.1) in cylindrical polar coordinates. In the present work, a central difference method was used with 551×201 grid points with a uniform mesh in *z* and a stretched mesh in *r*. The simplest approach is to neglect the eigenfunction structure and put the forcing into a boundary condition on *r*=*a*, where *a* is the radial location of the inflection point in the velocity profile. The source terms in (1.1) are set to zero within the solution domain; thus the flow is only forced by the imposed boundary condition. Following from considerations of the jump conditions for the Lilley–Goldstein equation (see Sandham *et al*. 2006*a*, appendix B), a pressure gradient boundary condition is appropriate and is applied at *r*=*a* by setting(2.8)where *ρ*_{0} is ambient fluid density. This equation is solved in a domain of size −25≤*z*≤40 and *a*≤*r*≤30. The simulation is run to time *t*≈30 and results are shown in a restricted subdomain −15≤*z*≤25 and *a*≤*r*≤15. This allows the initial transient to leave the subdomain, but no time for any reflections from the boundaries of the full domain to arrive back in the subdomain, thus external boundary effects are avoided. Unless otherwise noted, the velocity, density and sound speed profiles in the Lilley–Goldstein equation are taken at a fixed downstream location *z*_{0}=5.

## 3. Results

### (a) Velocity profiles and spatial growth factors

An analytic velocity profile was matched to the available experimental data of Stromberg *et al*. (1980). The profile gives the axial velocity at any *z* and *r* locations using(3.1)with empirical coefficients(3.2)and(3.3)Examples of the resulting velocity profiles are shown in figure 2, illustrating how the jet changes from a top-hat profile at *z*=0 to a fully developed flow downstream. Profiles of the radial component of velocity are computed using steady continuity equation and are used in the non-parallel PSE formulation. The jet is swirl free and so . For this base flow pressure is assumed to be uniform and temperature is obtained from a Crocco–Busemann relation(3.4)taking the jet temperature equal to the ambient temperature. The jet Reynolds number was set to *Re*=3600 and the jet Mach number to *M*=0.9 to match the conditions of Stromberg *et al*. (1980) and Freund (2001).

Note that the prescribed profiles are taken to fit experimental mean flow data that include the effect of turbulent stresses. Thus, one type of nonlinearity has already been included implicitly. In more general situations one could consider taking the base flow from a Reynolds-averaged Navier–Stokes calculation.

The PSE calculations continued downstream to *z*=20. The box size was chosen to be *L*_{r}=30 with 80 points in *r* and 81 points in *z*. From the output of PSE, a growth (*N*-) factor is defined by(3.5)Figure 3 shows *N*-factors computed using PSE for azimuthal modes *n*=0, 1, 2 and 3. In each case frequencies are chosen from *ω*=0.2 to 6.0 in steps of 0.2. From grid resolution studies, numerical errors in the peak *N*-factors are less than 1%. The peak *N*-factor is 5.9 for *n*=0 at *z*=4.75. The *n*=1 mode is nearly equally amplified (peak *N*-factor 5.6), while for the *n*=2 and 3 modes the peak amplification factor drops to 5.3 and 4.2, respectively. It is noteworthy that the *n*=1 mode gives the largest growth factors in the downstream (*z*>10) part of the jet. A jet Strouhal number is defined as *St*=*fd*/*U*_{J}, where *f*=*ω*/(2*π*) is the frequency, *d* is the jet diameter and *U*_{J} is the jet core velocity. For the *n*=0 mode the greatest amplification is found for *St*=0.45. This is in close agreement with Stromberg *et al*. (1980) who measured *St*=0.44 at *z*=4.2 and also commented that *n*=0 and 1 were approximately equal in the near field. At their maximum growth factor (near *z*=5) the higher azimuthal modes *n*=1, 2, 3 have *St*=0.35, 0.38, 0.38, respectively.

### (b) Source strength and basic radiation pattern

The sound source strength can be investigated by considering the amplitude of the mode interactions. Here we focus on the difference mode interactions that are found to be dominant at low radiation frequencies. To illustrate the procedure, figure 4*a* shows contours of |*A*^{−}| computed according to (2.5) for difference mode interactions of axisymmetric modes (*n*_{k}=0, *n*_{m}=0) for a range of frequencies *ω*_{j} (horizontal axis) and *ω*_{l} (vertical axis) at *z*=4.25. For *ω*_{j}=*ω*_{l} the corresponding wavenumbers are also equal, giving a small value for |*A*^{−}| (though not identically zero, since there is still a growth rate term in the pre-multiplying factor). The peak values of |*A*^{−}| are located at (*ω*_{j}=2.2, *ω*_{l}=3.4) and (*ω*_{j}=3.4, *ω*_{l}=2.2). The difference mode Δ*ω*=1.2 corresponds to *St*=0.19. Figure 4*b–d* shows similar contours for the (0,1), (1,1) and (0,2) interactions. The (1,1) interaction gives a similar structure compared with (0,0), and again has *St*=0.19. Asymmetric behaviour is seen for the (0,1) and (0,2) interactions. Larger values of |*A*^{−}| are consistently seen for a high-frequency axisymmetric mode interacting with a low-frequency azimuthal mode than for a low-frequency axisymmetric mode interacting with a high-frequency azimuthal mode. This is because the growth rates for axisymmetric modes reduce rapidly as the frequency is reduced, leading to lower growth factors for their interactions.

Figure 5*a* shows the maximum |*A*^{−}| over all frequency combinations, as a function of *z*. The (0,0) interaction peaks at *z*=4 and then decays rapidly. The (0,1) interaction has a higher peak value and a similar decay trajectory to (0,0), while the (1,1) interaction has a lower peak but remains level for increasing *z*, emerging with the highest |*A*^{−}| for *z*>8. It will be shown later that higher modes, for example (0,2), have low radiation effectiveness and so their amplitudes are not included here. The variation of the corresponding *St* with *z* is shown in figure 5*b*. Higher frequencies, up to *St*=1, are observed near the nozzle, but for *z*>2 the data level off, with 0.16<*St*<0.26. At *z*=4 (the location with the peak growth factor), we have *St*=0.19 for the (0,0) mode interaction. This can be compared to the far-field dominant frequency at 30° to the downstream direction (i.e. close to the peak radiation direction), found to be *St*=0.2 by Stromberg *et al*. (1980) and by Papamoschou (2007) for his single-stream jet case at *M*=0.9 and *Re*=3.2×10^{5}.

The structure of the boundary condition defined by equation (2.8) is shown in figure 6 for the (0,0) mode interaction with maximum |*A*^{−}|. Figure 6*a* gives the boundary condition as a function of *z* at *t*=0, while figure 6*b* shows the energy spectrum of this signal. It can be seen that, while the dominant peak is at a phase speed close to 0.5, there is a tail of the spectrum extending out to supersonic (*c*_{ph}>1.11) phase speeds. It is this part of the spectrum that is responsible for the sound radiation.

Figure 7*a–c* shows the radiation pattern arising from the mode interaction leading to the largest |*A*^{−}| for the (0,0), (0,1) and (1,1) mode interactions, respectively. The same contour levels have been used for each interaction. The (0,0) and (1,1) interactions both give axisymmetric radiation, predominantly in the downstream direction. By contrast the (0,1) mode gives only weak sound radiation. These results are in broad agreement with those of Stromberg *et al*. (1980), who measured the azimuthal variation of the phase of the spectral component with *St*=0.22 and found it to be predominantly axisymmetric. All interactions, and the (1,1) interaction in particular, show evidence in the near field for a response at the forcing frequency. This response is located downstream of the main wave packet input, typically for *z*≥12, and does not appear to radiate strongly into the far field. This response is not present if one sets in equation (1.1).

### (c) Radiation effectiveness

For each interacting mode pair, a measure of the effectiveness of conversion of the forcing at *r*=*R* into sound is given by(3.6)Here the numerator is proportional to the radiated sound power, with the integral taken over the angle *ψ*, measured relative to the downstream direction, at a sphere radius , the origin for which is taken at *z*=5 and *r*=0. To avoid near-field contributions, the integration range is taken from 30° to 150°. The denominator is proportional to the intensity of the input forcing, with the integral carried out over all *z*. An overbar indicates an average over one period of the wave motion and the quantity has been scaled, so that the effectiveness is unity for radiation from a spherical source in the high-frequency limit *ωR*/*c*_{0}→∞. Table 1 shows the numerator (in dB; note that only relative levels are important) for various mode interactions, while table 2 shows the radiation effectiveness *η*. These data are all given for the difference interaction of a mode with frequency *ω*_{j}=3.2 (together with *n*_{k} for each row) with a mode with frequency *ω*_{l}=1.8 (with *n*_{m} for each column), taking . This interaction (*St*=0.22), while not quite at the maximum for any of the difference mode interaction cases shown in figure 4, is close to the maximum for all interactions, with the (0,0) mode being the furthest from its maximum and hence underrepresented. In the far field the radiation is dominated by the (*n*_{k}=1, *n*_{m}=1) difference mode interaction. It can be seen that the radiation intensity and effectiveness of off-diagonal elements decrease rapidly as the difference azimuthal wavenumber increases. Figure 8 shows the directivity of the sound radiation ( at ), using the same frequency combinations as figure 7. It can be seen that the radiation at *ψ*=90° is 15–20 dB lower than that at *ψ*=30° and there is very little upstream radiation computed using this approach. In this case, the (0,0) interaction is comparable to (1,1), in contrast to the results in table 1, since it was taken at the frequency combination that gave the maximum |*A*^{−}|.

The sensitivity of the results to the choice of *z*_{0} is shown in figure 9, for *z*_{0} equal to 0, 5 and 10. As *z*_{0} is increased the shear layer thickness grows and the radiation in the forward arc is increased. For these cases, the lower limit of the computational domain for the wave equation is given from equation (3.2) for the particular choice of *z*_{0}. If, instead, *a* is kept fixed, then all the cases show the same radiation at *ψ*=90°, the value for which can also be found by solving a simplified problem with and variations in *ρ* and *c*_{0} neglected. Overall the variations seen in figure 9 amount to 3 dB at *ψ*=30° due to the choice of *z*_{0} in equation (1.1).

## 4. Discussion and conclusion

The simplified difference mode interaction model presented here is successful in predicting a component of subsonic jet noise arising from hydrodynamic instabilities, which may be efficiently computed with the PSE methodology. Some aspects of the results are in agreement with the experiments of Stromberg *et al*. (1980), notably the Strouhal numbers of the primary jet instability and of the radiated sound from the dominant interactions. The model shows that the axisymmetric and first azimuthal modes are active within the jet, while the sound radiation is primarily axisymmetric. The apparent origin for the sound is at *z*=5 (in units of jet diameter), which is a little upstream of the *z*≈7 seen in the visualizations of Freund (2001).

The low Reynolds number transitional jet considered here may of course be a pathological case where such a theory can work quite well, but not be applicable to higher Reynolds number applications. However, the experimental observations of instability waves within high Reynolds number jets (Suzuki & Colonius 2006) suggest that there is likely still to be a role for the sound produced by interactions between instability waves in fully turbulent jets. The current mechanism only applies to subsonic jets. At high Mach numbers the sound radiation comes directly from the primary instability (Tam & Morris 1980; also shown clearly in the PSE of Cheung *et al*. 2007). Salgado (2007) carried out linear PSE for a range of Mach numbers. Direct sound radiation in the downstream arc could be observed in pressure contour plots down to *M*=1.4, which may explain why Stromberg *et al*. (1980) observed a change of the dominant sound radiation frequency as the jet Mach number was changed from *M*=0.9 to 1.4. At higher Mach number, radiation Strouhal number coincided with near-field Strouhal number, in contrast to the low Mach number case where they were different.

It is worth noting that the present formulation allows the Lilley–Goldstein equation to be solved in the time domain. This is normally precluded by the development of instabilities within the shear layer. However, since we apply the forcing on *r*=*a* where *a* is the location of the inflection point in mean profile, there are now no inflection points within the domain and the problem of spurious unstable solutions is avoided.

The main limitations of the method given here are that the base flow has been taken to be fixed and the disturbances are assumed to be small. Thus there is no interaction between the base flow and the growing instabilities, which ought to be linked via Reynolds stresses. This effect can be brought into the PSE model, either using an energy equation as in Morris *et al*. (1990) or by performing a fully nonlinear PSE as in Cheung *et al*. (2007). It should be noted however that nonlinear PSE is not the same as the present nonlinear interaction approach. Nonlinear PSE still needs to be coupled to an acoustic analogy to produce the sound from a subsonic jet flow. Another consequence of the linearity of the present approach is that there is no phenomenon of vortex roll up, during which one frequency emerges as dominant and suppresses neighbouring frequencies. In the present model, mode interactions with small frequency differences lead to large sound radiation. This is despite the fact that the axial quadrupole amplitude *A*^{−} reduces near *ω*_{j}=*ω*_{l} (figure 4) and is due to the increased efficiency of the sound conversion for smaller frequencies. For example, a (0,0) mode interaction between *ω*_{j}=2.8 and *ω*_{l}=2.2 reached a conversion effectiveness of *η*=6% and radiates at 81 dB, well in excess of those shown in tables 1 and 2. In practice one would expect one or other of these frequencies to be suppressed as a vortex rolls up in the jet, with a result that such close interactions become less important. A possible simplification is only to consider pairing-type interactions, with factor of two differences between successive frequencies. Further studies using DNS of the Navier–Stokes equations are needed before one can say whether this effect can be included into the current approach. The increased radiation for smaller frequency differences also explains why the (0,0) mode in figures 6 and 7 had a higher radiation intensity than the other modes, whereas in tables 1 and 2 the (1,1) mode was dominant.

The present model can be extended (Sandham *et al*. 2008) to look at different source terms. Instead of introducing the forcing along the surface of a cylinder as done here, the various source terms on the r.h.s. of (1.1) can be computed from the PSE mode shapes and used in the numerical solution of a wave equation. The relative importance of the various sources can then be studied, with the potential for some insight into the physical mechanisms that lead to efficient sound production.

## Acknowledgments

This work is supported by EPSRC (grant no. EP/E032028/1) and Rolls-Royce plc. The authors would like to thank Prof. Chris Morfey, Dr Anurag Agarwal and Dr Zhiwei Hu for commenting on a draft version of the manuscript.

## Footnotes

One contribution of 10 to a Theme Issue ‘Theoretical fluid dynamics in the twenty-first century’.

- © 2008 The Royal Society