Solution of unsteady Euler equations: Gust–cascade interaction tones

Computers & Fluids 36 (2007) 724–741 www.elsevier.com/locate/compfluid

Solution of unsteady Euler equations: Gust–cascade interaction tones M. Nallasamy

a,*

, R. Hixon b, S. Sawyer

c

a QSS Group, Inc. NASA Glenn Research Center, Cleveland, OH 44135, USA Mechanical Engineering Department, University of Toledo, Toledo, OH 43606, USA c Mechanical Engineering Department, University of Akron, Akron, OH 44325, USA

b

Received 14 May 2005; received in revised form 19 March 2006; accepted 28 June 2006 Available online 10 October 2006

Abstract The problem of interaction of a vortical gust with a two-dimensional cascade is considered. Full nonlinear time dependent Euler equations governing the flow are solved employing a 6th-order accurate spatial differencing scheme and a 4th-order accurate time marching technique. The vortical gust is represented by a Fourier series which includes three harmonics. The acoustic response of the cascade for single and multi frequency (vortical) excitations are calculated. The solutions show the generation and propagation of modes that are expected from the theory. It is demonstrated that at low amplitudes of excitation, the time domain analysis produces characteristics of the propagating modes such as the complex mode amplitudes, phase variations, axial waveforms, and tangential waveforms that are in very good agreement with those expected from the linear theory. The exponential decay of the cutoff modes of the first harmonic is also clearly observed. The sound pressure levels of the propagating modes obtained from the present nonlinear time domain analysis are compared with the results of a linearized Navier–Stokes solution and a linearized Euler solution (frequency domain analyses) and good agreement between the results is observed for all the propagating modes. Ó 2006 Published by Elsevier Ltd.

1. Introduction The gust cascade interaction problem has been studied extensively using semi-analytical and numerical approaches due to its importance in understanding turbofan rotor stator interaction noise. Rotor stator interaction tone noise is the result of interaction of the periodic disturbances of the rotor blade mean wakes with the stator vanes. The rotor wakes impinging on the stator vanes produce fluctuating aerodynamic forces which generate noise. An accurate calculation of the fluctuating forces is the critical element in the interaction noise analysis and until recently such a calculation was done on a two dimensional basis. The stator row can be unrolled at a fixed radius, so that the stator becomes a linear cascade, thus reducing the problem to two dimensions and tractable. A two dimensional analysis is valid only when strong radial gradients *

Corresponding author. Present address: Mail Code ASRC-23, Kennedy Space Center, FL 32899, USA. E-mail address: [email protected] (M. Nallasamy). 0045-7930/$ - see front matter Ó 2006 Published by Elsevier Ltd. doi:10.1016/j.compfluid.2006.06.002

do not exit. Such an analysis is a good approximation for high hub-to-tip radius ratio (>0.75) turbofans [1]. Also, the three dimensional problem may be analyzed in terms of several two dimensional problems (so called strip theory [2,3]) and hence the extensive investigations on the two dimensional problem. The two dimensional gust cascade problem can be analyzed employing either a time linearized approach or a nonlinear time marching approach. A number of time linearized analyses have been carried out starting with the Sears problem of an airfoil encountering a gust. In a time linearized analysis, the gust is assumed to be a small harmonic perturbation of the uniform steady mean flow, yielding a system of equations in the frequency domain. Hence, in this approach typically the gust response to a single frequency excitation is studied [4–8]. In the linearized analysis, nonlinear effects cannot be assessed when the amplitude of excitation is no longer in the linear range. In the present time domain analysis, the full nonlinear time dependent Euler equations are solved numerically employing a 6th-order accurate spatial differencing scheme,

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

725

Nomenclature A interaction tone amplitude a1, a2, a3 excitation amplitudes C vane chord c1, c2, c3, c4 characteristic variables c speed of sound E, F flux vectors km mode eigenvalue Md duct Mach number m circumferential mode order Nb number of rotor blades Nv number of stator vanes n harmonic order Pi inflow stagnation pressure p static pressure po outflow static pressure Q conserved flow variable vector r radius r0 duct radius

and a 4th-order accurate time marching technique. The time domain analysis requires significantly larger computational time compared to the linearized Euler analysis. A time domain analysis has the advantage that the acoustic response to all the harmonics of interest can be extracted from one solution and thus may be able to mimic the real flow more closely. Also, linear/nonlinear regimes, self interaction, and interaction between frequencies [9] may be explored. In the nonlinear range, energy transfer between different frequencies (harmonics) occur, and such energy transfers are easily handled in a time domain approach. The computational time requirements and resolution constraints have so far limited the time domain analysis to single airfoil gust response [10] and flat plate cascade [11]. With the availability of parallel processing algorithms and computer clusters with computing power far excess of super computers, time domain approach has become feasible for gust cascade interaction study of real modern turbofan stage configurations. In this paper, the tones generated by the interaction of a vortical gust with a two dimensional cascade is studied using a time domain approach. The incident gust is the measured periodic mean rotor wake impinging on the stator vanes. The measured wake is represented by a Fourier series which includes only three harmonics of the blade passing frequency (BPF). A parallel computational aeroacoustic code developed for the gust–cascade interaction analysis is used in this study. The code employs a 6th-order spatial differencing scheme and a 4th-order time marching technique. Following the description of the gust cascade problem, the form of the Euler equations solved is presented. Brief descriptions and relevant equations employed for the spatial differencing scheme, time marching technique, and boundary conditions are given along

S Ti t V1 u, v ug, v g W1 x, y ai b c Dx Dt n, g k / x X

stator vane gap inflow stagnation temperature time absolute flow velocity mean velocity components in x, y directions gust velocity components in x, y directions relative flow velocity Cartesian coordinates inflow angle gust angle ratio of specific heats axial grid size time step curvilinear coordinates wavelength phase angle frequency rotational speed of the rotor

with a summary of the solution procedure. Interaction modes which are expected to be produced as a result of the gust cascade interaction are arrived at for the fan stage under consideration. The technique employed to extract modal information from the converged numerical solution is discussed and the computed acoustic response in the linear regime is presented in detail, discussed, and compared with the available linear analyses results.

2. Two dimensional cascade problem A sketch of the rotor stator interaction problem is given in Fig. 1a. The spinning/convecting rotor blade mean wakes impinge on the stator vanes generating what is called

Fig. 1a. Sketch of the rotor stator interaction problem.

726

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

the rotor stator interaction tone noise. The velocity triangle at the inflow to the stator describes the nature of the flow. The wake from the rotor blades (gust) is in the direction of relative flow velocity (W1) and is at an angle b to the xaxis. The direction of tangential velocity (Xr) is in the direction of rotation of the rotor as indicated, resulting in the absolute velocity (V1) making an angle a to the x-axis. Here X is the rotational speed of the rotor, and r is the fan blade radius. In the present study, the stator is modeled as an isolated linear cascade that is excited by the convected vortical gust of the rotor wake. Thus, the acoustic response of the stator is uncoupled from the rotor and can be independently computed. Fig. 1b shows the unrolled section of the stator at a fixed radius, to form a linear cascade of airfoils, giving a two dimensional representation of the rotor stator interaction problem studied in the present investigation. The cascade has a gap-to-chord ratio of S/C = 2/3, where C is the vane chord and S is the vane gap. The farfield boundaries (inflow and outflow planes) are located at 1.5C upstream and downstream from the center of the airfoil axial chord. The mean flow conditions at the inflow and outflow planes are given as

where Pi and Ti are the normalized mean stagnation pressure and stagnation temperature, respectively. ai is the mean flow angle and po is the mean static pressure at the outflow plane. The inflow periodic wake disturbance is described at the inflow plane as ~ u0g ðy; tÞ ¼ fa1 cosðk y y  xtÞ þ a2 cosð2ðk y y  xtÞÞ þ a3 cosð3ðk y y  xtÞÞg^eb ¼ 0;

p0g ðy; tÞ ¼ 0

^eb ¼ cosðbÞ^ex  sinðbÞ^ey ;

b ¼ 44

x ¼ 3p=4;

a1 ¼ 0:005;

a2 ¼ 0:003;

k y ¼ 11p=9; a3 ¼ 0:0007:

3. Governing equations The full nonlinear Euler equations governing the two dimensional cascade flow are solved in the present investigation. In Cartesian coordinates these equations are written as oQ oE oF þ þ ¼0 ot ox oy

ð2Þ

where 2

Inflow: Pi = 1, Ti = 1, and ai = 36° Outflow: po/Pi = 0.92

q0g ðy; tÞ

Here ^eb ¼ ðcos b;  sin bÞ is the unit vector in the direction of the relative flow. ^ex ; ^ey are the unit vectors in the x and y directions. x is the fundamental reduced frequency, ky is the transverse wavenumber, and ai’s are the gust harmonic amplitudes. This was a benchmark problem in the fourth computational aeroacoustics workshop [12]. The frequency is normalized by the chord divided by the ambient speed of sound, the wavenumber is normalized by the vane chord, and the gust amplitudes are normalized by the ambient speed of sound. The harmonic amplitudes of excitation a1 (at BPF), a2 (at 2BPF) and a3 (at 3BPF) have been varied to study the linear/nonlinear behavior of the generated tones.

ð1Þ

3 q 6 qu 7 6 7 Q¼6 7 4 qv 5 2

E qu

ð3Þ 3

6 qu2 þ p 7 6 7 E¼6 7 4 quv 5 uðE þ pÞ 2 3 qv 6 quv 7 6 7 F ¼6 2 7 4 qv þ p 5

ð4Þ

ð5Þ

vðE þ pÞ and   1 2 2 p ¼ ðc  1Þ E  qðu þ v Þ 2

ð6Þ

These equations are cast in generalized curvilinear coordinates employing the chain rule formulation as oQ on oE og oE on oF og oF þ þ þ þ ¼0 os ox on ox og oy on oy og

Fig. 1b. Two dimensional cascade problem representation.

ð7Þ

Eq. (7) is solved employing high order spatial differencing and time marching techniques. In the present investigation a multi-block non-overlapping grid is used. The grid metrics calculations require spatial derivatives of the grid variables. The same higher accuracy spatial differencing scheme (compact 6th-order, see below) is used to obtain the spatial derivatives of the grid variables. The grid generated using the Gridpro

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

(a commercial grid generation code) has provided reasonably smooth grid for the present problem. 4. Spatial differencing, time marching, and boundary conditions A brief description of the spatial difference scheme, time marching technique, and far field boundary conditions used in the current solution procedure is given below. Filtering scheme used to filter unresolved waves to avoid contamination of the solution is also discussed. 4.1. Spatial differencing A standard 6th-order difference form is written as       1 of  of  4 of  þ þ 5 ox iþ1 ox i1 5 ox i   1 1 7 ðfiþ2  fi2 Þ þ ðfiþ1  fi1 Þ ¼ Dx 60 15

i

on the gust–cascade interaction problem. The schemes comparison study showed the following: The explicit 2ndorder scheme produced wrong interaction tones and spurious modes while the explicit 6th-order scheme produced the expected modes along with spurious modes. The higher order schemes, optimized DRP and compact 6th-order were found to produce the interaction tones and their amplitudes accurately [14]. In this paper, the solutions were obtained with the compact 6th-order scheme. 4.2. Time marching method A 2N-Low Storage [15] optimized Runge–Kutta scheme [16] is used and it is called the 5–6 scheme because it alternates between five (m = 1) and six (m = 2) stages per time step. At each time step the flow variables are updated as follows:

ð8Þ

To solve for the derivative, a scalar tridiagonal matrix must be solved. To avoid this and to reduce stencil size from five points to three, a prefactorization is used which splits the above equation into forward and backward biased stencils which are solved separately and then added together ! F ! a of F of   þ 1  a ox iþ1 ox i   1 b 1b ðfiþ1  fi Þ þ ðfi  fi1 Þ ¼ Dx 2ð1  aÞ 2ð1  aÞ ! ! B a of B of   ð9Þ þ 1  a ox i1 ox i   1 b 1b ðfi  fi1 Þ þ ðfiþ1  fi Þ ¼ Dx 2ð1  aÞ 2ð1  aÞ  F B    of  of of ¼  þ  ox  ox ox i

727

i

where 1 1  pffiffiffi 2 2 5 a b ¼ a2  30 a¼

At the boundaries explicit seven point stencils are used to retain the 6th-order accuracy of the spatial differential scheme at the boundaries. Thus, the accuracy at the boundaries is 6th-order as it is in the interior. Forward and backward stencils are used at the start and end of compact difference sweeps [13]. At the block boundaries also the formal 6th-order accuracy of the spatial differencing scheme is maintained. This is done by employing an eleven point stencil and using data from the adjacent blocks. Four spatial differencing schemes, (i) explicit 2nd-order (ii) explicit 6th-order (iii) optimized DRP and (iv) prefactored compact 6th-order, available in the code, were tested

Ql ¼ al;m Ql þ Dt

oQ ot

Ql ¼ Ql þ bl;m Ql

ð10Þ

tl ¼ t0;m þ cl;m Dt where al,m, bl,m, and cl,m are constants. The method is 4thorder accurate and is stable to a CFL number of about 1.25. 4.3. Filtering Nonlinear equations are solved on arbitrary grids and some waves become unresolved as the solution proceeds. To prevent the unresolved waves from contaminating the solution and destroy its accuracy they are removed by filtering while taking care not to damp the resolved waves. The filter is used only on the conserved variables before calculating the fluxes. The interior filter stencil is written as 1 f i ¼ fi  ½252f i  210ðfiþ1 þ fi1 Þ 1024 þ 120ðfiþ2 þ fi2 Þ  45ðfiþ3 þ fi3 Þ þ 10ðfiþ4 þ fi4 Þ  ðfiþ5 þ fi5 Þ

ð11Þ

The scheme is 10th-order as described by Kennedy and Carpenter [17]. As we approach the boundaries, one-sided stencils are used and the reduced ordering of filter results. For the 10th-order filter employed in the code, an eleven point stencil is used. But at the boundaries and up to five points from the boundary, the filter is reduced to 5th-order accuracy. From the sixth point (and beyond) from the boundary, the filter is of 10th-order. 4.4. Nonreflecting boundary conditions Nonreflecting boundary conditions derived using the characteristics theory [18], are employed at the farfield boundaries. The characteristics are defined as

728

8 9 2 2 c1 > c > > > > > < c2 = 6 0 6 ¼6 > > c 3 > 40 > > > : ; c4 0

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

38 0 9 q> 0 1 > > > > = < u0 > qc 0 7 7 7 0 v > qc 0 1 5> > > > ; : 0> p qc 0 1 0 0

ðc1 Þt ¼ vðc1 Þy ð13aÞ

At the outflow, a 2nd-order formulation shown below is employed. ðc4 Þt ¼ uðc2 Þy  vðc4 Þy

ð13bÞ

A 4th-order formulation such as 1 ðc4 Þt ¼  ðc þ uÞðc2 Þy  vðc4 Þy 2

ð13cÞ

for the out flow boundary condition was found to produce numerical instability and hence only Eq. (13b) was used. The implementation of the nonreflecting boundary conditions is done as follows. The contribution to the characteristics due to the gust and the derivatives of the gust characteristics are written as ðcm þ um Þ ðcm  um Þ ðc3g Þy  ðc4g Þy 2 2 ðcm  um Þ ðc2g Þy  vm ðc3g Þy ðc3g Þt ¼  2 ðc4g Þt ¼ qm cm ððug Þt þ ðvg Þt Þ ðc2g Þt ¼ 

ð14Þ

The gust contributions to the characteristics are subtracted c2t ¼ ðc2 Þt  ðc2g Þt c3t ¼ ðc3 Þt  ðc3g Þt c4t ¼ ðc4 Þt  ðc4g Þt

  1 1 1 ðc ðc ðc Þ þ Þ þ Þt 1 t 3 t 4 c2 2 2

ðuÞt ¼

1 ½ðc3 Þt  ðc4 Þt  2qc

ð12Þ

The four characteristics c1, c2, c3, and c4 correspond to entropic, vortical, downstream propagating acoustic, and upstream propagating acoustic waves, respectively. The primed quantities are the perturbation variables. For a subsonic flow, the changes in the incoming characteristics are computed and outgoing characteristics are left unchanged. The incoming characteristics at the inflow boundary are c1, c2, and c3 while at the outflow boundary c4 is the only incoming characteristic. The inflow characteristics (time derivatives) are computed from a 4th-order formulation

1 1 ðc2 Þt ¼ vðc2 Þy  ðc  uÞðc3 Þy  ðc þ uÞðc4 Þy 2 2 1 ðc3 Þt ¼  ðc  uÞðc2 Þy  vðc3 Þy 2

ðqÞt ¼

ð15Þ

The derivatives of the flow variables are then computed as follows

ð16Þ

1 ðvÞt ¼ ½ðc2 Þt  qc 1 ðpÞt ¼ ½ðc3 Þt þ ðc4 Þt  2 At the inflow, gust components are now added ut ¼ ðuÞt þ ðug Þt vt ¼ ðvÞt þ ðvg Þt

ð17Þ

Once change in velocity components are computed, change in conserved variables, (qu)t, (qv)t, Et can be computed. 4.5. Mean flow boundary conditions To apply the mean flow boundary condition (MFBC), the mean flow is directly computed at the boundaries as the flow evolves. The time derivative of the flow at the boundary is written as a sum of three components: (Qt)boundary = (Qt)MFBC + (Qt)gust + (Qt)non-reflective. The method of evaluating (Qt)MFBC is discussed in detail in [19]. It is briefly stated here. For subsonic flows, the desired mean flow conditions are specified in terms of the inflow stagnation temperature and pressure, inflow angle and outflow static  t Þ is simply the flow variables pressure. The mean flow ðQ integrated over a given length of time. The mean flow integration is performed using the same integration routine in the code. At any point on the inflow boundary, the current mean flow error is first computed. Using 1-D characteristics, the incoming characteristics (c1, c2, c3) that will correct this error can be calculated as in [19]. The fourth characteristic is outgoing and it cannot be modified. Similarly, at the outflow, the first three characteristics are outgoing and cannot be modified. Only a change in c4, needed to correct the mean flow error, at the outflow is computed. Once the modifications to the characteristics are known, the necessary modification to the flow variable (DQ) can be computed. Then the first term, the mean flow boundary condition is  tÞ written as ðQ MFBC ¼ rðDQÞ. Here r is a relaxation factor used to speed convergence to the desired mean flow. In this work, a value of 1.0 was used for r, while efforts continue to find an optimum. The second term, (Qt)gust, specification is done as follows. We remove the incoming gust disturbances from the Giles condition (characteristics, (Eq. (15))) and then add the time derivative of the gust components in at the inflow boundary (Eq. (17)). The nonreflecting boundary conditions, (Qt)non-reflective, used in the code are the Giles inflow (Eq. (13a)) and outflow (Eq. (13b)) boundary conditions. This solution procedure has been found to be reliable and accurate to study the incident gust response.

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

729

5. Summary of solution procedure The full nonlinear time dependent Euler equations governing the two dimensional cascade flow are solved employing a time domain approach using the spatial differencing, time marching, and boundary condition implementations described above. These form the core of a parallel computational aeroacoustic (CAA) code [20] employed in the present study. Flow equations written in chain-rule curvilinear form are solved using a prefactored 6th-order compact scheme for spatial differencing. The time marching uses a 2N Storage 4th-order nonlinear extension of 5–6 Low Dissipation and Dispersion Runge–Kutta (LDDRK) scheme. An explicit 10th-order constant coefficient artificial dissipation is used at every stage of the Runge–Kutta solver to dissipate unresolved waves. On the airfoils, the time derivative of the velocity normal to the wall is set to zero. To apply the mean flow boundary condition (MFBC), the mean flow is directly computed at the boundaries as the flow evolves. At the inflow and outflow planes nonreflecting boundary conditions for the unsteady flow are implemented as described above. 6. Computational domain and numerical solution The fan stage considered in this paper has 22 rotor blades and 54 stator vanes (Fig. 1a). In the present solution procedure, because of symmetry, it is enough if we consider half the stator (circumference). That is, we need to solve the unsteady flow only for 27 vane passages. The computational domain is shown in Fig. 2a. A periodic boundary

Fig. 2b. Regular grid, 1.5 < x < 1.5: single passage.

condition in the y-direction is specified. The inflow outflow planes were first located at a distance of 1.5C from the origin as mentioned in the problem definition section. For the cascade configuration considered here, the number of interaction modes that are expected to be generated, the modes that propagate and the modes that decay among the generated modes are discussed in Section 7. The grid used initially in this study is a regular grid, 1.5 < x < 1.5 (Fig. 2b) consisting of eight grid blocks per passage with 9506 grid points in each passage. Fig. 2b shows four vane passages to indicate the cascade solidity. Care is taken to increase the grid resolution in the leading edge, trailing edge, and wake regions. The grid density has been found to be sufficient to produce resolutions required for the excitation frequencies considered in this study. With the specified gust amplitudes (a1, a2, and a3), the solution is run until periodicity in pressure is achieved on the airfoil and on the inflow and outflow boundaries. The solution was run up to 235 periods so that change in circumferential mode amplitudes between successive periods is within one percent. Once the periodic solution is established, the solution is processed to obtain the response of the acoustic modes generated. To avoid small reflections from the downstream boundary influencing the downstream propagating mode characteristics a stretched grid, 1.5 < x < 12.85 (Fig. 2c) was used for all the results presented in this paper as discussed in Section 9. 7. Gust–cascade interaction tones

Fig. 2a. Computational domain.

Excitation at the first three harmonics of the blade passing frequency is expected to produce response at these frequencies. The acoustic tones generated as a result of the impingement of rotor wakes on the stator vanes of a fan stage is given by the Tyler–Sofrin [21] criterion. Each

730

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

Fig. 2c. Stretched grid, 1.5 < x < 12.85: single passage.

harmonic tone may have one or more circumferential modes that propagate or decay. The idea of analyzing the fan tone noise in terms of spinning pressure patterns in the fan duct was first introduced by Tyler and Sofrin. If the number of rotor blades is Nb and the number of stator vanes is Nv then the circumferential modes generated is given as Circumferential mode order m ¼ nN b  kN v where, n is the harmonic order and k is an integer. The circumferential mode m represents the number of waves that can fit in the circumference and thus an integer. The modes generated at the stator leading edge propagate in the upstream (x) and downstream (+x) directions. A circumferential mode order m at a frequency x will propagate only if its cutoff ratio is greater than unity. If the cutoff ratio is less than unity the mode is called a ‘‘cutoff’’ mode and decays exponentially away from the stator leading edge in the upstream and away from the stator trailing edge in the downstream directions. The cutoff ratio, n of a mode m is defined as n¼

ðx=cÞr0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : k m 1  M 2d

ð18Þ

Here km is the mode transverse eigenvalue, c is the speed of sound, r0 is the fan duct radius, and Md is the duct axial Mach number. The fan stage considered in this paper has 22 rotor blades and 54 stator vanes (Fig. 1a). The modes generated as a result of the rotor wake stator interaction at the expected response frequencies is given by Tyler–Sofrin criterion. The modes generated and their cutoff ratios are shown in Table 1. A positive sign indicates that the mode

is spinning in the direction of rotation of the rotor, while a negative sign indicates spinning opposite to that of rotor rotation. The cutoff ratio as defined in (18) has km which depends on the hub to tip radius ratio. Here the cutoff ratio is calculated using the upstream duct Mach number, and the hub to tip radius ratio near the stator leading edge as in [2,22]. In the present problem, the hub to tip radius ratio is constant in the upstream duct. At the first harmonic (n = 1), the generated modes have cutoff ratios less than one (Table 1). Thus we have the opportunity to examine the exponential decay of the cutoff modes at the first harmonic. At the second and third harmonics, the modes generated have cutoff ratios greater than one and we can study their propagation characteristics. However, we note that the cutoff ratio of mode 42 is very close to unity and this gives rise to interesting propagation characteristics as will be discussed in the results section. Because of symmetry, it is enough if we consider half the stator (circumference), (i.e.) Nv = 27 vanes (passages) corresponding to Nb = 11 as mentioned in the section above. Accordingly, the mode numbers will be half as shown in Table 1; they are at 1BPF 11 and 16, at 2BPF 5, and at 3BPF 6 and 21.

8. Joint temporal–spatial transform Once the solution is converged, the cascade unsteady flow as a function of time is known over a specified period. The acoustic pressure response is not only harmonic in time, but also periodic in the tangential direction y. From the known acoustic pressure p as a function of (x, y, t), we need to extract the acoustic mode amplitude P as a function of axial location, spatial mode order, and

Table 1 Circumferential mode orders and cutoff ratios Harmonic order, n

Nb = 22, Nv = 54 mode, m

Nb = 11, Nv = 27 mode, m

Cutoff ratio, n

Cutoff/propagates

1 1 2 3 3

22 32 10 12 42

11 16 5 6 21

0.64 0.454 2.673 3.399 1.051

Cutoff Cutoff Propagates Propagates Propagates

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

frequency (x, m, n). This is done by applying a joint temporal–spatial transform N 1 NX Dt 1 2 X P ðx; m; nÞ ¼ pðx; y l ; tj Þe2pinj=N Dt e2piml=N NN Dt l¼0 j¼0

ð19Þ where, Dt is the time between samples, NDt is the number of samples over interval T, frequency f = n/T, and n = 0, 1, 2, . . . (NDt  1)/2 is the harmonic number. N is the number of points in the tangential direction. m is the spatial (circumferential) mode order, N/2 < m < (N/2  1). Note that only integer values of m are permissible, corresponding to the number of periods that fit in the circumference and including only positive frequencies implies that the rotor and its convected gust are rotating in the positive direction. To perform the Fast Fourier transform (FFT), the number of time samples employed is usually a power of two. Since FFT is done over the entire computational domain, the number of samples cannot be too large. But, it was found that a minimum of 32 samples (point in the time series) is needed to avoid aliasing problems. No discernible change in amplitude was observed when the sample size was increased to 64.

731

9. Results and discussions The cascade acoustic response for a vortical disturbance is studied using the time domain approach. The acoustic response for the multi frequency excitation (at the first three harmonics) with amplitudes, a1 = 0.005(BPF), a2 = 0.003(2BPF), and a3 = 0.0007(3BPF), is described. Linearity of the response for a single frequency excitation at high amplitudes is explored. 9.1. Steady airfoil pressures and harmonic loading The steady airfoil loading obtained from the current time marching procedure is compared with that obtained from a nonlinear Euler solver called TURBO [23] in Fig. 3a. The agreement between the two solutions is excellent and only small discrepancies are observed at the trailing edge region which may stem from the way trailing edge point is treated in these two solution procedures. The oscillations near the leading edge in Fig. 3a and b stem mainly from insufficient airfoil surface definition near the leading edge for the grid generator and this was corrected in the subsequent grid generations. Next we look at the magnitude of the unsteady pressures. For the low amplitudes of

Fig. 3. Airfoil surface pressures: (a) steady surface pressure comparisons with Turbo code solutions; (b) unsteady pressure variation on the airfoil surface, minimum and maximum variations; (c) airfoil surface pressure variations at first three harmonics.

732

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

excitation considered here, the unsteady (total) loading itself is a small percentage of the steady loading (Fig. 3b) as would be expected for a linear problem. The loading at each harmonic is shown in Fig. 3c. The loading decreases with increase in harmonic order for this cascade. 9.2. Propagating modes at 2BPF and 3BPF: excitation at three harmonics For the cascade–gust interaction problem considered here, the theory predicts propagation of m = 5 mode at 2BPF and two modes m = 6 and m = 21 at 3BPF. However, m = 21 mode is very close to cutoff (Table 1). From the converged solution of known acoustic pressure p as a function of (x, y, t), the acoustic mode amplitude P as a function of axial distance, mode order, and harmonic order is obtained as described above (in Joint temporal–spatial transform). The amplitude variations along the axial direction of the propagating modes are plotted in Fig. 4 for the inflow and outflow regions. (The extent of the inflow region 1.5 < x < 0.8, and the outflow region 0.8 < x < 1.5 are indicated in Fig. 2b). In Fig. 4, two sets of results are presented: one from the regular grid (1.5 < x < 1.5, Fig. 2b) solution and the second from the stretched grid (1.5 < x < 12.85, Fig. 2c) solution. First we focus on the regular grid solution (dashed lines in Fig. 4). We see that the amplitude of mode m = 5 at 2BPF is constant along x in the inflow region as expected of a propagating mode. However, in the outflow region, propagating wave m = 5

shows signs of a small reflection from the downstream boundary (The mode amplitude deviates slightly from the constant value). An acoustic wave of spatial mode order m at harmonic n traveling downstream from the stator is expected to behave like: P ðx; m; nÞ ¼ P dn ðm; nÞeikx;dn x where P dn ðm; nÞ is the amplitude of the downstream propagating wave. In the presence of an upstream going wave reflected from the outflow boundary, the transformed pressure would be modulated such that P ðx; m; nÞ ¼ P dn ðm; nÞeikx;dn x þ P up ðm; nÞeikx;up x where P up ðm; nÞ is the amplitude of the upstream going wave. The reflection was found to be approximately three percent of the incident wave amplitude for this mode with the regular grid solution. The acoustic response at 3BPF shows that m = 6 mode propagates clearly in the upstream and downstream directions (Fig. 4). Mode m = 21 propagates only in the upstream direction and decays slowly in the downstream direction. A cutoff mode would decay exponentially along x (as in Fig. 12 for the modes in BPF). The proximity of m = 21 to cutoff and small reflections from the downstream boundary seemed to contribute to the observed behavior of the mode m = 21 in the downstream (Fig. 4b). To avoid those small reflections, it was decided to stretch the grid in the downstream direction to x = 12.85 and obtain the solution. Before going to the stretched grid solution, a solution on a fine grid having nearly three times the number of grid points (27, 294 grid points per passage (not shown)) but with the same axial extent (1.5 < x < 1.5) was obtained. A comparison of

Fig. 4. Amplitudes of the propagating modes: (a) upstream propagating modes, inflow region; (b) downstream propagating modes, outflow region; (c) airfoil pressure distribution at 3BPF: comparison of regular and stretched grid solution.

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

the mode sound pressure levels (SPL in dB computed as 20 log (prms/pref), where pref = 20 lPa), from the regular

733

and fine grid solutions, at the inflow and outflow boundaries showed good agreement (Table 2) except for the

Table 2 Grid study: sound pressure levels (SPL): inflow and outflow planes Mode order, m

Inflow (x = 1.5) SPL (dB)

Outflow (x = 1.5) SPL (dB)

Coarse

Fine

Change

Coarse

Fine

Change

1BPF

16 11

105.46 96.74

105.68 95.93

0.22 0.81

98.00 113.24

100.0 113.37

2.00 0.07

2BPF

5

122.74

122.57

0.14

123.75

123.65

0.10

3BPF

6 21

106.25 103.23

105.94 103.28

0.31 0.05

102.12 103.68

101.69 102.93

0.43 0.75

(Coarse grid: 9506 points; fine grid: 27,294 points).

Fig. 5. Features of acoustic response at second harmonics: (a) temporal transform, real part of unsteady pressure; (b) spatial transform, circumferential mode orders; (c) circumferential mode waveforms in tangential direction.

734

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

decaying mode m = 16 at the outflow. The comparison showed that: (i) a fine grid (with 27,294 grid points per passage) did not affect the behavior of the downstream going m = 21 mode, and (ii) the amplitudes of the propagating modes did not show appreciable change as shown in Table 2, (iii) the regular grid with 9506 points per passage had 124 grid points in the axial direction providing enough resolution (points per wavelength) for the frequencies considered and thus produced accurate results. The

minimum number of grid points per wavelength is 30 at the highest frequency of excitation, 3BPF, which is enough to produce the required resolutions. Hence, in the stretched grid in Fig. 2c, the grid density in the region 1.5 < x < 1.5 is retained as in Fig. 2b. In the rest of the study, to avoid the small reflections from the downstream boundary influencing the propagation characteristics of the downstream going waves, a grid stretched gradually in the downstream direction to

Fig. 6. Features of acoustic response at third harmonics: (a) temporal transform, real part of unsteady pressure; (b) spatial transform, circumferential mode orders; (c) circumferential mode waveform in tangential direction at inflow; (d) circumferential mode waveform in tangential direction at outflow.

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

735

Fig. 7. Complex pressure amplitudes of propagating modes: (a) propagating modes at 2BPF; (b) upstream propagating modes at 3BPF, inflow region; (c) downstream propagating modes at 3BPF, outflow region.

x = 12.85 (Fig. 2c) was employed. The stretched grid is intended to dissipate the outgoing waves at large distances in the downstream direction. The solution obtained with this grid was analyzed and the mode amplitudes are included in Fig. 4 (solid lines). First, we see that the mode m = 21 also clearly propagates (indicated by the constant mode amplitude along x) in the downstream direction (Fig. 4b), as expected from the linear theory. The propagating mode amplitudes from the regular and stretched grid solutions for the other modes compare well at the inflow and outflow regions. At the inflow, the two solutions produce nearly similar amplitudes of the three propagating modes. At the outflow, m = 5 and m = 6 mode amplitudes are the same for the two solutions. The effect of small reflections on the axial variation of 2BPF amplitude in the outflow region (dashed line) is not seen with the stretched grid solution and mode m = 21 clearly propagates in the downstream direction. Also, note that in the downstream direction, the mode amplitudes remain constant much beyond the original downstream boundary location of x = 1.5. No spurious modes (interaction tones) were revealed in the spatial–temporal transforms. We now examine the stretched grid solution in detail.

The pressure distributions on the airfoil surface obtained with regular and stretched grids for the 3rd harmonic are shown in Fig. 4c. The plots show substantial differences (on this expanded scale) in pressure along the airfoil chord between the two solutions. While the

Fig. 8. Axial waveforms at the inflow for the propagating modes, m = 5, 6, and 21 at 2nd and 3rd harmonics.

736

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

differences in magnitude are seen through the entire chord, the characteristics are different only near the leading edge region. Even here, the leading edge peaks are nearly identical. But beyond the leading edge peak pressure differences occur. In this region, the pressure level on the pressure surface is higher for the stretched grid solution while for the regular grid solution the suction surface pressure level is higher (Fig. 3c). It appears that this change in loading near the leading edge region is instrumental in changing the downstream propagation characteristics of m = 21 mode. Since the solutions obtained with the stretched grid is found to give mode characteristics that are in excellent agreement with those of propagating modes of the linear theory, the stretched grid has been employed for all the solutions examined in this paper. 9.2.1. Tangential waveforms and spatial–temporal transforms Spatial, temporal transforms and tangential waveforms of the second harmonic response are shown in Fig. 5. The temporal transform, (the real part of acoustic pressure), shown in part (a) of the figure. clearly shows that the mode m = 5 propagates in the upstream and downstream directions. (Note that the solution for 27 passages is shown). The spatial transform (part b) shows the amplitude variation along the axial direction of the only propagating mode m = 5. The tangential waveform for this mode may be written as   2p p5 ¼ A5 cos y ð20Þ k5 where A5 is the mode amplitude and k5 is the tangential wavelength. We see that the waveforms extracted from Fig. 5a for the inflow and outflow regions plotted in part Ó of the figure shows the expected waveforms. Since we

Fig. 9. Phase variations of mode pressures along the axial direction at the inflow for propagating modes, m = 5, 6, and 21 at 2nd and 3rd harmonics.

have only one propagating mode, the waveforms for the inflow and outflow regions are similar and differ only by the difference in mode amplitudes and phase difference. The acoustic response at the 3rd harmonic has two propagating modes m = 6 and 21. The spatial–temporal transforms and tangential wave forms produce interesting features shown in Fig. 6. The interactions of the two propagating modes produce different mode patterns at the inflow and outflow. We observe from Fig. 4 that at 3BPF, mode m = 6 dominates in the inflow region while at the outflow region mode m = 21 dominates. Accordingly, the temporal transform in Fig. 6a, shows the mode 6 at the inflow region and the mode m = 21 at the outflow region. The spatial transform shows two modes propagating in the upstream and downstream directions (Fig. 6b).

Fig. 10. Axial waveforms at the outflow for the propagating modes, m = 5, 6, and 21 at 2nd and 3rd harmonics.

Table 3 Axial wavelength n

m

Inflow region, kx

Outflow region, kx

2 3 3

5 6 21

0.8 0.63 0.60

1.90 1.30 7.82

Fig. 11. Phase variations of mode pressures along the axial direction at the outflow for propagating modes, m = 5, 6, and 21 at 2nd and 3rd harmonics.

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

The tangential waveforms in the inflow and outflow regions may be expressed as     2p 2p ð21Þ pin ¼ A6in cos y þ A21in cos y þ /in ; k6 k 21    2p 2p pout ¼ A6out cos y þ A21out cos y þ /out ; ð22Þ k6 k21 respectively. Fig. 6c shows the tangential waveform extracted from the present nonlinear time marching solution (Fig. 6a) in the inflow region compared with the

737

analytical waveform of Eq. (21). The inflow mode amplitudes of A6in = 0.00005 and A21in = 0.0002 (Fig. 4a) were used in the equation. A value of /in = 0.75p was used to align the waveforms. The agreement of the waveform from the time marching solution with the analytical waveform of Eq. (21) is very good. The tangential waveforms at the outflow are shown in Fig. 6d. Here the waveform from the time marching solution is compared with the analytical waveform of Eq. (22). The outflow mode amplitudes of A6out = 0.00003 and A21out = 0.000092 (Fig. 4b) were

Fig. 12. Cutoff mode (BPF) characteristics: (a) inflow region; (b) outflow region; (c) variation of complex pressure amplitudes along the axial direction.

Fig. 13. Sound pressure levels, (SPL) in dB, of the Acoustic modes generated as a result of gust interaction with the cascade at x = 1.5 (inflow) and x = 1.5 (outflow).

738

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

used. A value of /out = p/4 was used to align the waveforms. The agreement between the two waveforms is again found to be extremely good. The inflow and outflow region tangential waveforms clearly bring out the nature of the interacting modes, m = 6 and 21, resulting in the observed waveforms. 9.2.2. Axial waveforms, complex amplitudes, and phase variations of propagating modes Fig. 7 shows the complex amplitude of the propagating modes in the inflow and outflow regions. A perfectly propagating wave would appear as a circle in such a plot. At 2BPF (Fig. 7a), m = 5 mode propagates well in the upstream and downstream directions. The mode propagating in the downstream direction has higher amplitude than that of the upstream propagating one. In Fig. 7b, complex amplitudes of the two modes (m = 6 and 21) at 3BPF are shown for the inflow region while Fig. 7c shows for the outflow region. The perfect circles in the two regions clearly indicate well propagating modes. The mode m = 6 dominates in the inflow region while the mode m = 21 dominates in the outflow region as shown in Figs. 4 and 6. A propagating acoustic mode generated as a result of rotor stator interaction is the combination of a spinning mode and a traveling axial wave. The circumferential waveforms of the propagating modes were presented above. In

this section, we present the axial waveforms and the mode phase variations. The acoustic modes generated at the stator leading edge propagates in the upstream (x) and downstream (+x) directions. The acoustic mode traveling in the upstream direction travels against the flow direction and thus the axial wavelength is reduced. On the other hand, the mode traveling in the downstream direction travels along the flow direction and thus the axial wavelength is increased. This behavior is illustrated here: The axial wave forms of the propagating modes in the upstream direction are shown in Fig. 8. Linear phase distributions of the upstream propagating modes (in the inflow region) of m = 5, 6, and 21, are shown in Fig. 9. The axial waveforms of the modes propagating downstream (outflow region) are shown in Fig. 10, which clearly show the increased wavelengths (Table 3) compared to the wavelengths of the upstream propagating waves in the inflow region. Fig. 11 shows the phase distribution of downstream traveling waves (in the outflow direction) of modes m = 5, 6, and 21. A linear phase variation is clearly observed for all the modes as would be expected from the linear theory. Thus, the present nonlinear time marching solution is shown to produce clearly the expected characteristics of all the upstream and downstream propagating modes as illustrated in amplitude, phase, axial waveform, and tangential waveform plots.

Fig. 14. Comparisons of the sound pressure levels in dB of the propagating modes in the upstream and downstream directions with two linearized solutions: (a) 2BPF: m = 5; (b) 3BPF: m = 6; (c) 3BPF: m = 21.

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

9.3. Decaying (cutoff) modes at BPF: excitation at three harmonics For the cascade problem considered here, the blade passing frequency is cutoff (Table 1). The BPF modes that are generated at the stator vane leading edge (m = 11, called the rotor locked mode and m = 16, the interaction mode) are expected to decay exponentially in the upstream

739

and downstream directions. This behavior is illustrated in Fig. 12a for the inflow region and in Fig. 12b for the outflow region, where the regular and stretched grid solutions show exponential decay of these modes. A look at the complex amplitude variation of the modes (Fig. 12c) show random variations of the amplitudes as opposed to perfect circles of the propagating modes (Fig. 7).

Fig. 15. Linear behavior of single frequency excitation at BPF: (a) inflow, m = 11; (b) outflow, m = 11; (c) inflow, m = 16; (d) outflow, m = 16.

Fig. 16. Linear behavior of single frequency excitation at 2BPF, m = 5: (a) inflow region; (b) outflow region.

740

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

9.4. Sound pressure levels The sound pressure levels of the modes generated at the three harmonics at upstream (x = 1.5) and downstream (x = +1.5) locations are shown in Fig. 13. The mode m = 5 of the second harmonic has the maximum sound pressure level consistent with the experimental measurements [24] and the level at the outflow is about 1.1 dB higher than that at the inflow. A direct comparison of the absolute sound pressure levels cannot be made since the excitation frequency and amplitude were modified to ‘better’ define the benchmark problem. Also, the experimental data incorporate three dimensional effects and hence direct comparisons will be made when the 3-D solutions (currently in progress) become available. Of the two propagating modes in 3BPF, the mode m = 6 dominates in the inflow while the mode 21 dominates in the outflow. In the cutoff BPF, the mode 16 has a higher sound pressure level than that of mode 11 at the inflow, while the reverse is true at the outflow. Note that at these measurement locations (inflow and outflow boundaries) the cutoff modes have amplitudes comparable to those of the 3rd harmonic propagating modes. The sound pressure levels of the propagating modes, m = 5, 6, and 21, at the upstream (x = 1.5) and downstream (x = +1.5) locations are compared with those of two frequency domain analyses in Fig. 14. The results labeled Linflux [23] were obtained using a frequency domain linearized Euler code called LINFLUX which has been tested at NASA Glenn Research Center and United Technology Research Center on rotor-stator inter-

action problems. The results labeled Hydlin [25] were obtained using a linearized Navier–Stokes code, HYDLIN. This code has been tested for a range of turbomachinery unsteady problems including noise generation and propagation at Rolls Royce. The sound pressure levels obtained from the present time marching approach agrees with those of Hydlin for all the propagating modes. The Linflux results for the 3rd harmonic mode m = 6 differ substantially from the present results as wall as those of Hydlin, due to the grid resolution problems. A fine grid LINFLUX solution for 3BPF is not available at this time for comparison. 9.5. Single frequency excitation: effect of amplitude of excitation In this section, acoustic response characteristics of the cascade for single frequency excitations at BPF, 2BPF, and 3BPF at high amplitudes of excitation are explored. The acoustic response amplitude for the single frequency excitation at BPF, a1 = 0.005 (Fig. 15) is the same as found above for the three-harmonic excitation. Next, the amplitude of excitation is increased 10 times to a1 = 0.05 and the response computed for this excitation amplitude is also shown for the inflow and outflow regions in Fig. 15. In the figure, the dashed lines represent the response for the excitation amplitude a1 = 0.005 multiplied by a factor of 10. The linearity of the response and decay characteristics of the modes are clearly observed except for the mode m = 16 in the outflow region which shows deviations from the linear behavior.

Fig. 17. Linear behavior of single frequency excitation at 3BPF: (a) inflow, m = 6; (b) outflow, m = 6; (c) inflow, m = 21; (d) outflow, m = 21.

M. Nallasamy et al. / Computers & Fluids 36 (2007) 724–741

Fig. 16 shows amplitudes of the propagating mode m = 5 for the 2BPF excitation in the upstream and downstream directions. When the excitation amplitude is increased 10 times to a2 = 0.03, the 2BPF response also increases 10 times as shown in the figure. The dashed lines represent the response for the excitation amplitude a2 = 0.003 multiplied by a factor of 10. The linear response of the propagating mode for the 2BPF excitation is clearly exhibited. For the 3BPF excitation, a similar increase in excitation amplitude by a factor of 10 to a3 = 0.007 (from the original 0.0007) is shown to increase amplitudes of the propagating modes m = 6 (Fig. 17a and b) and m = 21 (Fig. 17c and d), by the same amount in the inflow as well as in the outflow regions. Again, the dashed lines represent the response for the excitation amplitude of a3 = 0.0007 multiplied by a factor of 10. The acoustic response of the cascade for the 3BPF excitation at a3 = 0.007 is clearly seen to be linear. High frequency excitation in this work is limited to ten times the original amplitude. The upper limit up to which the response remains linear has not been established. However, for single frequency excitation at the blade passing frequency, at ten times the original amplitude, nonlinear effects begin to appear (Fig. 15d) indicating that we exceeded the linear limit. This high amplitude of excitation at blade passing frequency was found to modify the response at higher harmonics. The nonlinear effects at single frequency excitations will be reported in a future paper. Multi frequency excitations at high amplitudes produce interaction between frequencies and self-interaction [9] leading to nonlinear response of the generated modes. The time domain analysis provides a unique opportunity for studying the nonlinear interactions. The results of the high amplitude multi frequency excitations will be presented in a future paper. 10. Conclusions A time domain analysis was carried out to examine the acoustic response of a two dimensional cascade for an incident vortical disturbance. Full nonlinear time dependent Euler equations governing the flow are solved using high order accurate spatial differencing and time marching techniques on parallel processors. We show that at low amplitudes of excitation, the characteristics of the propagating modes obtained from the time domain analysis are in very good agreement with those expected from the linear theory. The computed sound pressure levels agree with those of the linearized Navier–Stokes (Hydlin) solution. The results of single frequency excitation show linear acoustic response over a wide range of excitation amplitudes. Work is in progress to extend the time domain analysis to compute rotor stator (three dimensional) interaction tones.

741

References [1] Morfey CL. Rotating pressure patterns in ducts: their generation and transmission. J Sound Vib 1964;1:60–87. [2] Envia E, Nallasamy M. Design selection and analysis of a swept and leaned stator concept. J Sound Vib 1999;228:793–836. [3] Nallasamy M, Envia E. Computation of rotor wake turbulence noise. J Sound Vib 2005;282:649–78. [4] Smith SN. Discrete frequency sound generation in axial flow turbomachines. Aeronautical Research Council report and Memoranda No. 3709; 1973. [5] Whitehead DS. A finite element solution of unsteady two dimensional flow in cascades. Int J Numer Meth Fl 1990;10:13–24. [6] Hall KC, Verdon JM. Gust response analysis for cascades operating in nonuniform mean flows. AIAA J 1991;29:1463–71. [7] Hall KC, Clark W. Linearized Euler predictions of unsteady aerodynamic loads in cascades. AIAA J 1993;31:540–50. [8] Verdon JM. Linearized unsteady aerodynamic analysis of the acoustic response to wake/blade row interaction. NASA/CR-2001210713; 2001. [9] Morse PM, Ingard KU. Theoretical acoustics. New York: McGrawhill; 1968. [10] Lockard DP, Morris PJ. Radiated noise from airfoils in realistic mean flow. AIAA J 1998;36:907–14. [11] Hixon R, Shih SH, Mankbadi RR. Evaluation of boundary conditions for gust cascade problem. J Propul Power 2000;16:72–8. [12] Proceedings of the fourth computational aeroacoustics workshop on benchmark problems. NASA-CP-2004-212954; 2004. [13] Hixon R. Nonlinear comparison of high-order and optimized finite difference schemes. J Comput Fluid Dyn 2000;13:257–77. [14] Hixon R, Nallasamy M, Sawyer DS, Dyson R. Comparison of numerical schemes for a realistic computational Aeroacoustic benchmark problem. Int J Aeroacoustics 2005;4:379–97. [15] Stanescu D, Habashi WG. 2N-Storage low dissipation and dispersion Runge–Kutta scheme for computational acoustics. J Comput Phys 1998;143:674–81. [16] Hu FQ, Hussaini MY, Manthey JL. Low-dissipation and lowdispersion Runge–Kutta scheme for computational acoustics. J Comp Phys 1996;124:177–91. [17] Kennedy CA, Carpenter MH. Several new numerical methods for compressible shear layer simulations. Appl Num Math 1994;4: 397–433. [18] Giles M. Nonreflecting boundary conditions for Euler equation calculations. AIAA J 1990;28:2050–8. [19] Hixon R, et al. Mean flow boundary conditions for computational aeroacoustics. AIAA Paper 2003-3299; 2003. [20] Hixon R, Nallasamy M, Sawyer SD. Parallelization strategy for computational aeroacoustic code. AIAA paper 2002-2583; 2003. [21] Tyler LM, Sofrin TG. Axial flow compressor noise studies. Trans SAE 1962;70:309–32. [22] Nallasamy M. Computation of noise radiation from fan inlet and aft ducts. J Aircraft 1997;34:387–93. [23] Envia E. Bench mark solution for the Category3-problem2: cascade– gust interaction. In: Fourth computational aeroacoustics on benchmark problems. NASA/CP-2004-212954; 2004. p. 59–65. [24] Woodward RP, Hughes CE, Jeracki RJ, Miller CJ. Fan noise source diagnostic test – Farfield acoustic results. AIAA paper 2002-2427; 2002. [25] Coupland J. Linear unsteady CFD analysis of the cascade–gust interaction problem. In: Fourth computational aeroacoustics workshop on benchmark problems. NASA/CP-2004-212954; 2004. p. 97– 105.