Measurements and computational fluid dynamics predictions of the acoustic impedance of orifices Journal of Sound and Vibration

The response of orifices to incident acoustic waves, which is important for many engineering applications, is investigated with an approach combining both experimental measurements and numerical simulations. This paper presents experimental data on acoustic impedance of orifices, which is subsequently used for validation of a numerical technique developed for the purpose of predicting the acoustic response of a range of geometries with moderate computational cost. Measurements are conducted for orifices with length to diameter ratios, L/D , of 0.5, 5 and 10. The experimental data is obtained for a range of frequencies using a configuration in which a mean (or bias) flow passes from a duct through the test orifices before issuing into a plenum. Acoustic waves are provided by a sound generator on the upstream side of the orifices. Computational fluid dynamics (CFD) calculations of the same configuration have also been performed. These have been undertaken using an unsteady Reynolds averaged Navier – Stokes (URANS) approach with a pressure based compressible formulation with appropriate characteristic based bound- ary conditions to simulate the correct acoustic behaviour at the boundaries. The CFD predictions are in very good agreement with the experimental data, predicting the correct trend with both frequency and orifice L/D in a way not seen with analytical models. The CFD was also able to successfully predict a negative resistance, and hence a reflection coefficient greater than unity for the L = D ¼ 0 : 5 case.

engine). If the amplitude of the incident wave is not high enough to induce reverse flow through the liner, the absorption is typically linear (i.e. the absorbed energy relative to the incident acoustic energy is constant). In such a case the ratio between the fluctuations of pressure and flow rate across the perforated liner is a constant complex number. This is generally referred to as the acoustic impedance, and in the linear regime, this is independent of the amplitude of pressure fluctuation. This paper is concerned with the linear absorptive properties of a range of aperture geometries with bias flow in the presence of plane incident waves.
An approach is required that accurately captures the flow mechanisms and absorption properties of aperture geometries that are typically used within engineering applications. For the approach to provide a useful design tool it must be capable of providing results within a relatively short time frame (e.g. days) and, in the long term, must be capable of applying to a whole range of aperture shapes and sizes. In this paper experimental and computational methods are presented that have been used to obtain the acoustic impedance of circular orifices subjected to incident plane acoustic waves over a range of frequencies. The results from these methods are then presented for comparison. This is done for apertures of length to diameter ratio 0.5, 5 and 10. This is thought to reflect the range of geometries likely to be encountered in practice. It includes cases in which the flow forms a vena-contracta downstream of the aperture (L=D ¼ 0:5) to flow reattachment inside the aperture, with the flow then remaining attached over a significant portion of the aperture length (L=D ¼ 10:0). The numerical method employed here uses unsteady Reynolds Averaged Navier-Stokes (URANS) modelling to produce time resolved simulations of the unsteady flow induced by the acoustic waves. Before describing the experimental and numerical methods in detail, we first discuss the relative advantages of the numerical methods available.
It is typically assumed that the apertures within a porous plate are sufficiently separated so they behave as if operating in isolation. In this way the various investigations undertaken on a single orifice can be used to help understand the performance of perforated liners. For example Howe [2] developed an expression for the Rayleigh conductivity of an aperture through which a high Reynolds number bias flow passed. The theoretical model provided an exact analytical solution for an infinitely thin aperture. This was based on the harmonic pressure difference across the aperture causing the unsteady shedding of vorticity from the aperture rim. In this way acoustic energy is converted into turbulent flow energy that is subsequently dissipated as heat. This has been used as the basis for many predictions relating to perforated liners, for example [3,4] and Rupp and Carrotte (2011) [5]. However, for practical engineering applications the aperture must be of finite thickness (or length) and so an additional term was typically included. As indicated by [3] the thickness of the aperture (or its effective acoustic length) can be crucial in determining the performance of the aperture. However, it should be acknowledged that in many engineering applications the aperture length will significantly affect the flow field that is inherent to the model developed by Howe. For example, downstream of an infinitely thin orifice the flow will form a venacontracta, but at greater orifice length to diameter ratios (e.g. L=D 42) the flow reattaches inside the aperture. The mechanisms by which acoustic energy is absorbed will therefore be affected and the addition of a single term to capture this may be somewhat simplistic. Alternative strategies to those based on Howe include, for example, the use of a momentum balance across the aperture such as that outlined by Bellucci et al. [6]. In this case the model requires an aperture loss (or discharge) coefficient, along with an effective length, to enable the prediction of damping performance. The limitations associated with both of the above approaches are reflected in the fact that in conjunction with an effective length, the discharge coefficient that provides the best match with experimental acoustic absorption measurements does not often agree with the experimentally measured discharge coefficient obtained from the mean flow [7]. Zhong and Zhao devised a time-domain numerical tool for acoustic damping of a perforated liner with bias flow by transforming the frequencydomain analytical models into the discrete time-domain through z-transform [8]. Satisfactory agreements with the experimental data were obtained. However, as this method is based upon a 1D acoustic wave model, a suitable analytical model must already be known. This approach will have the limitations associated with such a model in that its use is restricted only to those flows for which the assumptions behind the model are valid.
CFD methods offer an attractive solution for predicting acoustic behaviour, as CFD calculations of mean flow fields are already a routine part of engineering design. Linearised methods for calculating acoustic behaviour have been used, for example Gikadi et al. [9] used this method to investigate a swirl burner. They used a technique whereby a RANS simulation was produced of the steady field before a method using Linearised Navier-Stokes Equations (LNSE) was used to find the magnitude and phase of reflected and transmitted waves through a swirler. The LNSE part of the calculation was carried out using a finite element method using a much coarser mesh than the RANS calculation.
Time resolved CFD methods offer the ability to reveal details of the unsteady flow as well as removing some modelling assumptions. Orifice flows of engineering interest will typically be at Reynolds numbers where the flow is turbulent. This turbulence will both determine the behaviour of the mean flow and play a key role in dissipating the energy put into the flow by acoustic forcing. As with the simulation of any turbulent flow a choice must be made as to how to treat this turbulence. At one extreme is the approach of using Direct Numerical Simulation (DNS) to resolve all temporal and spatial scales of turbulence. Zhang et al. performed DNS to study the acoustically excited flows through a circular orifice backed by a hexagonal cavity [10]. Accurate results were obtained up to a frequency of 2.5 kHz. Tam et al. also used DNS in their numerical investigation of resonators in three dimensions [11]. By virtue of its high-fidelity non-modelling approach, DNS is able to reproduce the small scale details of acoustic phenomena (such as dissipation mechanisms) accurately, even for high frequency ranges. However, the computational resources required by DNS limited both Zhang and Tam's works to only a small number of data points. For studies in lower frequency ranges (e.g. the one considered in the present work) the computational cost is expected to be excessively high. Another popular choice for academic research is Large Eddy Simulation (LES) in which the smaller spatial scales of turbulence are modelled rather than resolved. Andreini et al. developed a compressible Large Eddy Simulation (LES) technique [12] which has been used to study the acoustic damping properties of the multi-perforated plates used as cooling devices in gas turbine combustors. The CFD method also allowed investigation of the effect on the acoustic damping of geometrical details such as the staggering and angle of the holes [13]. Following the work of [4], Mendez et al. applied LES to analyse the damping effect of a perforated liner with bias flow [14]. In another recent study Lacombe et al. [15] also used compressible LES to study whistling behaviour of a single orifice from which they were able to achieve good agreement with experiment. Even though LES has been shown to give accurate results at a reduced cost compared to DNS, the simulations conducted in [15] still required weeks of calculation on 256 processors for a relatively simple geometry. The suitability of LES as a practical design tool, particularly for complex geometries, is thus significantly compromised.
A recent development in CFD simulations is the use of Lattice-Boltzmann solution methods, which offer advantages in terms of computational efficiency. These methods use a temporally and spatially resolved approach to simulating turbulent flows with the addition of a model to simulate the effect of spatially unresolved scales of turbulence, and so can be thought of as similar conceptually to the LES method. Lattice-Boltzmann methods have recently been used by da Silva et al. [16] to investigate the sound reflection at the open end of a duct issuing a subsonic mean flow, while Ji et al. [17] investigated excited flows through a small circular orifice. These works were able to find results which compared well with traditional LES techniques in terms of accuracy at a lower computational cost. These two studies were carried out for flows with Reynolds numbers of less than 1000, which is well below those typically found in industrial applications as investigated in this work.
URANS techniques use the approach of solving for the statistical average of the turbulent flow using some model to account for the effect of the turbulence on this averaged flow. Studies of whistling orifices using incompressible RANS techniques with a suitable post-processing technique have been performed [18,19] and report that the main physical features of the whistling phenomenon are reasonably reproduced. The incompressible nature of this method places restrictions on Mach number and the geometry size relative to the wavelength. These restrictions can be removed by using a compressible RANS formulation, although a density based compressible formulation will have the restriction that a CFL based on the speed of sound will limit the timestep. Pressure-based compressible methods offer the advantage that a larger flow speed based time-step can be used. Gunasekaran and McGuirk [20] used a URANS approach and developed a mildly compressible pressure based formulation and employed suitable boundary conditions to investigate the absorption coefficient of a duct-fed orifice. Recent work [21] developed a non-reflecting boundary condition for use with the fully compressible pressure based solver available in the open source OpenFOAM CFD code [22]. This was successfully used to predict the behaviour of a silencer. In this work we also use a fully compressible pressure based solver in OpenFOAM and, following the work in [20], employ suitable boundary conditions to simulate a siren delivering a specified signal at the inlet of the duct and a plenum on the other side of the orifice. In doing so we predict the variation with frequency of specific impedance for three aspect ratios of orifice using a relatively inexpensive computational technique applicable to more complex geometries. Whilst this work is concerned with wave propagation parallel to the hole centreline, it is thought that the methodology is suitable for propagation parallel to the hole but this would require further validation (see for example [4]). In addition, the proposed approach is also thought amenable to cavity backed resonant liners as studied in [14] where a relatively expensive LES formulation is used.
URANS methods offer a less expensive alternative for high Reynolds number flow as they can use coarser meshes and, crucially for acoustic work, they do not require long simulated times to average out the turbulent part of the flow. The disadvantage of using RANS techniques is that, by averaging the turbulence, the full process by which the energy is dissipated by the turbulence is not reproduced in the simulation. However if the RANS approach is able to successfully model the 'net result' of the turbulence on an acoustically forced flow field then its relatively low cost (compared to DNS or LES) has the potential to make it useful as an engineering tool to find the acoustic response of geometries where analytical models are not appropriate. To this end, its ability to simulate the impedance of orifices of L/D ¼ 0.5, 5, 10 is assessed by comparison to experimental data.

Experimental method
The experimental investigation was undertaken within an isothermal test facility operating at atmospheric conditions. The facility, shown in Fig. 1, is described in more detail by Rupp [7]. The basic test rig consisted of a 120 mm Â 120 mm square crosssectioned duct which was terminated at each end by large plenum chambers (50 m 3 ). At one end of the duct plates of varying thickness could be attached. Each plate could also incorporate orifices of various diameters. Hence, in this way a range of orifice L/D ratios could be investigated. Note that the number of orifices in each plate was also varied, this being a compromise between maximising the impedance being measured whilst ensuring effects due to radiation were maintained at an acceptable level. However, in all cases the orifice pitch was never less than 3.3 orifice diameters (i.e. so that each orifice could be assumed to act in isolation from its neighbour). A centrifugal fan attached to one plenum provided a mean pressure drop so as to pass air along the rectangular duct which exhausted to the downstream plenum via the orifice plate. For the results reported the mean pressure drop across the orifice plate (Δp=p) was approximately 0.5 percent (i.e. 500 Pa). This was measured via a static tapping located upstream of the orifice plate referenced to the pressure in the downstream plenum. Attached to the duct were 2 loudspeakers that generated plane acoustic waves which were varied in frequency from 50 Hz to 1 kHz (in 20 Hz increments). At the maximum frequency (1 kHz) it is still valid to assume plane axial waves within the duct, since the cut-on of higher order modes will only occur at frequencies in excess of approximately 1.4 kHz. All tests presented within this paper were undertaken within the linear absorption range with the acoustic pressure amplitudes within the duct being set to approximately 130 dB. The downstream exhaust plenum provides a steady (i.e. p 0 ¼ 0:0) exit boundary condition. Hence the unsteady (0 to peak) pressure drop across each orifice plate was approximately 20 percent of the mean pressure drop. However, linearity of the response was confirmed by tests being undertaken over a range of unsteady pressure amplitudes which showed that the acoustic characteristics (impedance, absorption etc.) were constant when made suitably non-dimensional. The combination of operating conditions and orifice geometries enabled investigation of orifice characteristics over a range of Strouhal numbers, defined as 2πfD=U, from 0.04 to 3.2 and L/D ratios from 0.1 to 10.

Multi-microphone technique
Measurements of unsteady pressure within the tube were obtained using 4 Kulite dynamic pressure miniature transducers. The transducers were statically calibrated by applying a known (i.e. reference) static pressure. Thereafter dynamic calibration was demonstrated by placing the transducers in the duct at the same axial location and, with the loudspeakers activated, all 4 transducers were shown to measure the same amplitude and phase over a range of frequencies. The transducers were connected to a National Instruments LabView Data acquisition system as described by [23,7]. The output from each transducer was suitably conditioned and acquired via a suitable 16 bit data acquisition card. The unsteady pressure signals were acquired at a frequency of 40 kHz in eight blocks of 32,768 samples (0.8192 s per data block). Each time history was Fourier transformed using the Fast Fourier Transformation method with a Hamming window. Thereafter an ensemble average over the 8 time histories was conducted on the Fourier transformed amplitudes. The auto-power spectra of each signal were calculated using the averaged Fourier transformed pressure amplitudes [24]. Furthermore one of the transducer signals was chosen as a reference to which all the other signals were compared. A cross-power spectrum then enabled the phase and magnitude of each unsteady pressure signal, relative to the reference value, to be determined.
At any point in the ductp the acoustic pressure is a superposition of the incident pressure wave travelling towards the plate (i.e. with the bias flow in this case) and the reflected wave away from the plate (upstream in this case). Ifp represents the complex amplitude of the waves theñ p ¼pexpðjωtÞ ¼p i expðjωt À jk þ xÞþp r expðjωt þ jk À xÞ where k 7 ¼ ω=ðc7UÞ and the subscripts 'i' and 'r' denote the incident and reflected waves respectively. The acoustic velocity can be found asũ If the pressure within the duct is measured at two locations then By determining the complex amplitude of the superimposed wave at multiple points, the incident and reflected wave magnitudes can be determined from the above expressions using the so-called two microphone technique [25]. The numerical predictions and experimental measurements were undertaken well below the cut-off frequency of the duct. In addition, the incident and reflected wave components were determined by interrogation of the flow field at sufficient distance from each orifice to ensure only the presence of plane waves i.e. any orifice local near field effects had decayed. In addition, the accuracy of this approach is strongly dependent on the axial location of the microphones relative to the pressure nodes and anti-nodes within the duct. Hence over the frequency range of interest this can lead to large errors and certain frequencies. In this case a multi-microphone method was therefore used [26]. Four dynamic pressure transducers were used along with the linear equation system Ax ¼b: This can be solved using a least squares fit whilst minimising the square of the residual J r With the complex amplitudes of the incident and reflected waves known, the acoustic pressure and velocity can be found at any position in the duct using Eqs. (1) and (2). Using a two microphone method it was found that pressure and wave amplitudes were in error by up to 3.3 percent and phase angles by up to 0.51. For the operating conditions being considered here this approximates to acoustic loss and absorption coefficients being accurate to within approximately 7 percent and 1 percent respectively. It is assumed that these errors are a worst case scenario and will be significantly reduced by application of the multi-microphone analysis method.

Orifice impedance
The acoustic characteristics of orifices can be described in a variety of ways. In this work we chose to use impedance to quantify the response. Impedance for the orifices can be defined as the complex ratio of acoustic pressure across the orifice to the unsteady part of the velocity induced through it: where R represents the in-phase component, known as resistance, and X is the out-of-phase component, known as reactance.
The unsteady pressure across the orifice can be found by assuming that the unsteady pressure at some point on the plenum side is zero while that on the duct side can be found from Eq. (1) with the x position set to that of the front of the plate. The unsteady velocity is found by evaluating Eq. (2) at the plate location and multiplying the amplitude by the area ratio of the duct to the orifice to account for the flow expansion. The location of the pressure node on the plenum node will not, in fact, be located at the interface between the plate and the plenum as the flow in the plenum in the vicinity of orifice will also be disturbed. This extra mass of air excited by the pressure wave is usually referred to as the 'end correction' and it is common practice to include the end correction in the inertia of the orifice, as is the case with Bellucci's model [6]. We have followed this practice by including the end correction in the calculated impedance; the unsteady pressure drop used in Eq. (6) is, in effect, that between the duct side and some point after the end correction.
The impedance based on the duct side unsteady pressure and the unsteady orifice velocity includes both the orifice and radiation impedance (i.e. the latter being associated with the sound radiated from the orifice). However, as described by Ingard and Issing [27] and Rupp [7] the reactive part of the radiation impedance is included in the orifice impedance. This is because it represents the effects of the inertial mass of the local air motion in the vicinity of the orifice (and hence is included via the end correction). The radiation resistance is similar to a plane piston in a rigid baffle. For the configurations considered here this is much smaller than the viscous resistance and so can be neglected. Hence the orifice impedance can be determined from the duct unsteady pressure (acting on the face of the orifice) and the orifice velocity. It should be noted that both experimental and numerical impedances have been calculated in the same manner and are, as such, directly comparable.

Analytical models for impedance
Analytical models are reported in the literature from which the impedance of an orifice can potentially be estimated. For example, one can obtain the impedance from the Rayleigh conductivity defined as [28] where A is the area of the orifice. In Howe's [2] well-accepted linear model the Rayleigh conductivity is determined as where I 1 and K 1 are the modified Bessel functions of the first and second kinds respectively and U v is the convective velocity associated with the shed vorticity (and usually assumed to be equal to the mean velocity through the orifice). Note that the Howe model is only applicable to infinitely thin orifices. However, other authors such as Jing and Sun [29] adapted it to orifices of finite lengths and expressed the specific impedance in terms of K R as where L is the length of the orifice. This leads to the addition of an extra term jρ 0 Lω to account for the orifice length. Hence for the case of zero length Eq. (9) reverts to what would be obtained from the original Howe model. This expression was applied and validated for length-to-diameter ratios up to 0.625 in [29]. However, for even higher L/D ratios such as those considered here (e.g. L/D¼5 and 10), the orifice is long enough for the separated bias flow to reattach to the orifice wall. This is likely to result in a substantial change to the orifice impedance (as will be subsequently discussed). For this reason, in the data presented here the Jing and Sun model is only applied to the orifice of L=D ¼ 0:5 and compared with the associated experimental and CFD datasets. An alternative approach was developed by Bellucci et al. [6] based on the momentum balance of a one-dimensional incompressible flow across the orifice. The semi-empirical model incorporates the discharge coefficient, in addition to the effective length, as an input parameter. If these parameters are known, the model is potentially applicable for both short and longer orifices. Therefore it is employed here for all the orifices considered. In a one-dimensional model, the relation between the velocity and pressure fluctuations can be expressed in general as whereũ andp are the complex amplitudes of velocity and pressure fluctuations respectively, L is the length of the orifice, L e is the end correction, ζ ac is the loss coefficient corresponding to the pressure loss due to area change at the orifice exit and ζ vis is the loss coefficient accounting for the viscous effect. When the incoming pressure wave is not strong enough to induce reverse flow through the orifice, ζ ac is given in terms of the discharge coefficient C d in Bellucci's derivation as where U c is the velocity at vena contracta. The discharge coefficient in Eq. (11) is obtained from experiment by passing a metered mass flow through the orifice and measuring the pressure drop across it. According to [30] the end correction of an orifice with an unlimited flange is 4D=3π and this is applied to both ends of the orifice. This leads to a total correction of 8D=3π [6]. The viscous loss, ζ vis , is expressed as [30] ζ vis ¼ where ν is the kinematic viscosity. The specific impedance is hence determined from Eq. (10) as The real part of the above equation represents the resistance while the imaginary part is the reactance. As the ffiffiffiffiffiffiffiffiffi ffi 2νω p will be small, this model will give a resistance that is nearly constant with frequency and a reactance that will increase linearly. Note also that use of this model requires a priori knowledge of the discharge coefficient of the orifice.

Computational method
A finite-volume URANS CFD method has been employed using the open-source arbitrarily unstructured CFD code OpenFOAM [22]. This uses the PISO algorithm which at each timestep iteratively solves the energy, momentum and compressible pressure correction equations in order to satisfy mass, momentum and energy conservation. Some details of the PISO algorithm are included below. In order to correctly reproduce the physical test case it is important that acoustic waves are correctly treated at the boundaries of the CFD domain. The advantage of implementing these into the OpenFOAM code is that other desirable practical features such as unstructured meshing and post-processing are already available. These will become more desirable if the method is applied to more complex geometry. The disadvantage of using an unstructured code is that higher order transport schemes such as the WENO scheme [31] are not available. However this disadvantage can always be offset with increased mesh resolution. For the frequencies of interest for thermoacoustic problems (less than 1 kHz when scaled to room temperature) the number of cells per wavelength suggested by [20] to give acceptable levels of dissipation with a second-order TVD scheme [32] are far from prohibitive in this work.

PISO algorithm
The Pressure Implicit with Splitting of Operator (PISO) algorithm is an intrinsically unsteady scheme and solves for the flow field at each timestep through a predictor step and a series of corrector steps. It has been shown to be suitable for flows across a range of speeds from low Mach number to supersonic [33]. The method is not new to this work and for a full description of the compressible formulation of the algorithm the reader is referred to [34], but we include some details here for clarity. At each timestep an initial predictor step is used in which conservation equations for momentum and internal energy are solved using the pressure field from the previous step. The updated internal energy is then used to update the thermodynamic state of the system including the compressibility defined for an ideal gas as ∂ρ=∂p ¼ 1=RT. This compressibility is then used in the subsequent iterative corrector loop. This loop consists of solving a transport equation, derived from mass conservation, for the pressure (rather than pressure correction as in the SIMPLE algorithm) using the current velocity field followed by correcting the velocity and density using the updated pressure. The corrector loop is run until convergence; i.e. the changes in flow properties have fallen below some desired levels at each iteration. The compressible nature of the flow is included through the use of the compressibility in the pressure equation. The use of a pressure based compressible scheme of this type allows for larger timesteps to be used than would be the case for a density based scheme where the CFL number would be based on the speed of sound rather than the flow velocity.

Boundary conditions
To give correct behaviour of acoustic waves at the domain boundaries it is necessary to use appropriate boundary conditions. In this work we have used the approach of local one-dimensional inviscid characteristic boundary conditions [35] that Gunasekaran and McGuirk [20] successfully applied to find the absorption coefficients of duct-fed orifices. The derivation of these boundary conditions starts from the 1D equations for the propagation of an acoustic perturbation in primitive variables derived by Hirsch [36]. For downstream and upstream travelling waves respectively these are where w down and w up are respectively the downstream and upstream travelling characteristics defined as If the boundary of the CFD domain is placed in a position where plane waves can be expected to cross the boundary then these 1D characteristics will be valid at the boundaries. By adding and subtracting Eqs. (14) and (15) and employing a firstorder Euler temporal discretisation scheme it is possible to derive equations for the change in fluctuating velocity and pressure in terms of the upstream and downstream characteristics: If w down and w up are known at the boundary then the fluctuating pressure and velocity at the boundary can be found. At a subsonic outlet, using a first-order upwind scheme, it is possible to find the outgoing characteristic w down using where the subscript 'i' indicates a face value whereas 'int' represents the first interior node. Similarly at a subsonic inlet the outgoing upstream characteristic is The remaining, incoming, characteristics can either be specified as a function of the outgoing characteristic or from a known incoming signal. This method allows for the treatment of a range of plane wave behaviours at the CFD boundary. By setting the magnitude of the incoming characteristic equal to that of the outgoing one a fully reflective boundary can be achieved. If w n in ¼ Àw n out then the reflection will be made to occur at a velocity node (i.e. u 0 ¼ 0). If w n in ¼ w n out then the reflection is made to occur at a pressure node (i.e. p 0 ¼ 0), as would be the case at the opening of a pipe to the atmosphere. The characteristic boundary condition can also be used to give the correct combination of a wave propagating out of the domain and another propagating into the domain. The incoming characteristic is based on the existing boundary values and also of the values of p 0 1 and u 0 1 from the incoming driving signal. For a downstream propagating sine wave these are p 0 1 ¼ A sin ð2πftÞ and u 0 1 ¼ p 0 1 =ðρ m c m Þ so that the incoming, downstream, characteristic at an inlet is and for an upstream propagating sine wave at an outlet they are p 0 1 ¼ A sin ð2πftÞ and u 0 1 ¼ Àp 0 1 =ðρ m c m Þ so that the incoming, upstream, characteristic at an outlet boundary is In this work the flow is essentially axial and with no non-acoustic variation with time at both inlet and exit of the CFD domain, the plenum exit boundary is more than 50D downstream of the orifice so any vorticity can be expected to have dissipated. As such the acoustic components calculated using Eqs. (20)- (23) can be added to the steady boundary conditions which give the desired mean flow field. Boundary conditions for pressure and velocity have been written in OpenFOAM to implement this. They have been written such that the fluctuations are superimposed onto the existing mean values of boundary pressure and velocity found in the calculation of the steady-state flow. The steady-state flow is first found by setting total pressure at the inlet and static pressure at outlet to give the desired pressure drop. The mean velocity at inlet and outlet will be found as part of the steady-state solution given these boundaries. The acoustic fluctuations are then added to these steady values at each timestep. The next subsection discusses how an outlet to a plenum can be simulated using characteristic based boundary conditions.

Plenum boundary condition
In the previous work of [20] using characteristic boundary conditions to investigate orifice absorption coefficients, an experimental set-up was used in which a duct was located both up and downstream of the orifice plate. In order to assume that the fluctuating pressure drop across the orifice is equal to the fluctuating pressure at the speaker side of the orifice, it is necessary to have a pressure node (p 0 ¼ 0) at the opposite side of the orifice to the noise generator. To do this with a duct fed arrangement, as in [20], it is necessary that the duct length must be matched to half the wavelength of the signal wave as there will be a pressure node at the open end of the duct. This is a major drawback when the aim is to investigate the impedance at a range of frequencies. Therefore the experimental set-up in this work is modified such that there is a plenum on the other side from the noise generator. As a result the exit of the orifice is exposed to a practically unlimited space and the pressure fluctuation downstream of the orifice can be assumed to be zero for the range of frequencies measured. This set-up presents a challenge to the CFD method to accurately simulate these conditions. The method employed here is based on recognising two properties of a plenum; firstly that the plenum represents a very large increase in flow area from the orifice such that any wave has to move a much larger mass of air when moving from orifice to plenum and so the amplitude is greatly reduced. Secondly any acoustic wave that does propagate into the plenum will not be reflected back. The non-reflecting behaviour can be reproduced using a characteristic boundary at the downstream outlet of the plenum where the incoming, upstream, acoustic signal is set to zero (i.e. in Eq. (23) p 0 1 ¼ u 0 1 ¼ 0) which gives the purely outward wave propagation with no reflection. To represent the effect of the expansion of the wave into the plenum it was found that if the CFD domain downstream of the orifice plate was made wider, the amplitude of the wave propagating downstream (i.e. through the plenum) was reduced. Beyond a certain point the acoustic energy lost into the plenum was small ($ 1 percent), and making the duct wider still beyond this point gave only a small change to the results obtained. Hence this width of 'plenum' was deemed suitable in this case.

Simulation details
The computational domain is illustrated in Fig. 2. Each of the plates tested in the experiments consists of nine orifices in order to improve the signal-to-noise ratio. The numerical simulation is simplified by considering one of the nine orifices only, i.e. the centre one specifically. Therefore upstream of the orifice a 40 mm Â 40 mm section of the duct corresponding to the part of the duct around the centre orifice is included in the domain, this is bounded by symmetry planes. The duct extends about five times its width upstream of the furthest pressure tapping in the experiment in order to ensure the quasi one-dimensionality of the acoustic field at all the measurement points. The plenum is modelled as a chamber with square cross-section with the side boundaries of this also set to symmetry planes. As discussed in Section 4.2.1 the plenum side of the CFD domain needs to be sufficiently large in order for only a small amount of acoustic energy to be transmitted upstream and give a pressure node at the upstream face of the orifice. A 120 mm Â 120 mm cross-section was used. Results for L=D ¼ 0:5 were repeated using a plenum of 240 mm with negligible change to results hence the geometry shown in Fig. 2 is deemed sufficient for an accurate reproduction of the actual scenario. The surface of the orifice plate is set to be a non-slip wall.
A multi-block 'O-grid' topology is used for the computational mesh. The centre blocks of the O-grids are formed with hybrid prismatic/hexahedral cells by sweeping the surface mesh through the domain from the inlet to outlet. The surrounding blocks are composed of hexahedral cells. A cross-section of the duct just upstream of the orifice is plotted in Fig. 3. On the non-slip walls, the first internal grid point is placed at yþ % 2 based on the flat-plate approximation. The maximum cell size in the duct is bounded in such a way that for the frequency range of concern, one wavelength is resolved by at least 50 cells. Meanwhile, the cell size at the outlet and boundaries of the plenum is allowed to grow further to provide additional numerical attenuation on the radiated waves. The total number of grid cells used ranges from approximately 1:0 Â 10 6 for L=D ¼ 0:5 to 1:5 Â 10 6 for L=D ¼ 10. As the acoustic waves have to be resolved the mesh resolution in the duct has to be relatively high and for the L=D ¼ 10 case, for example, only 20 percent of the mesh cells are in the orifice itself. This means that if the method were to be applied to a more complex geometry than a circular orifice the percentage increase of the total mesh size would not necessarily increase dramatically. A few indicators are used to examine the quality of the meshes. The equiangular skewness is always less than 0.5 and the lowest Jacobian ratio is 0.36. The average aspect ratio of the cells is about 2.8. A small number of cells with aspect ratio between 6 and 30 are found near the inlet and outlet of domain. These cells are deemed acceptable since the flow in these regions is mostly axial.
The mean bias flow is calculated prior to the unsteady acoustic flow. To calculate the mean bias flow, the total pressure and a velocity direction normal to the inlet plane are specified. The outlet boundary conditions are fixed pressure and advective velocity. A steady flow field is established using these boundary conditions before the acoustic simulation is commenced. The mean pressure drop was set to 500 Pa to correspond with the experiment. The transport properties are those standard for air at room temperature; i.e. ρ ¼ 1:22 kg/m 3 and ν ¼ 1:5 Â 10 À 5 m 2 /s. Using these the Mach and Reynolds numbers based on peak velocities in the orifice are 0.1 and 24,000 respectively. On account of its well-proven performance for free-shear and wall-bounded flows, the Menter k-ω shear-stress-transport (SST) model is adopted [37] for both steady and unsteady flows. This model allows the boundary layer to be resolved down to the viscous layer (i.e. low yþ). A secondorder TVD scheme [32] is used for spatial discretion in the momentum equations and turbulence equations, whilst firstorder upwind scheme is applied in the energy equation. For temporal discretisation the more stable first-order Euler backward scheme is selected. A PISO algorithm is used to solve the momentum and compressible pressure correction equations. The mean bias flow is calculated by running an unsteady solver until flow properties are seen to stop changing. A time step of 0.0005 s is used for calculation of the mean bias flow, mainly out of consideration of stability.
The acoustic simulations are run using the mean bias flow as the initial flow field. The acoustic boundary conditions described in the previous section are implemented at both the inlet and outlet. The acoustic signal is injected from the inlet using downstream travelling wave components p 0 1 ¼ A sin ð2πftÞ and u 0 1 ¼ p 0 1 =ðρ m c m Þ in Eq. (22). The amplitude A¼50 Pa and f is set to the desired frequency. The time step needs to be significantly decreased in the acoustic simulations in order to  reduce numerical dissipation on the acoustic waves. When the second-order central spatial scheme and first-order Euler backward temporal scheme are applied to the one-dimensional scalar wave equation, the attenuation of the wave over one wavelength can be obtained from modal analysis [38]. The time step in the acoustic simulations is so varied that the attenuation factor estimated in this way is comparable for the different frequencies of concern. It ranges from 1 Â 10 À 5 to 1:25 Â 10 À 6 s as the frequency increases. A full calculation of one geometry, including mean flow calculation and a sweep of ten frequencies, took approximately 1500 CPU hours which is relatively cheap by modern engineering computing standards. This makes the method attractive as a practical design tool.

Processing of results
Previous work [20] calculated absorption coefficients at a single frequency, showing the ability of the method to predict correctly the magnitude of the reflected pressure wave. In this work we use the CFD method to predict impedance to demonstrate its ability to predict both the magnitude and the phase of the acoustic response. In order to provide direct comparison between experiment and simulation impedances have been calculated from the CFD data in the same way as in the experiment (i.e. from analysis of the waves in the upstream duct section). Time series of pressure and velocity are taken from the CFD simulation at the position of the closest pressure tapping to the orifice plate. It was also confirmed in the simulations that at the tapping locations there was a plane wave, with the same signal being seen at the boundary and on the orifice centre line. These time series are taken over a time equal to 12 periods at the current frequency once a steady signal is achieved. Asũ andp are known from the time series, the up and down stream travelling pressure waves can be calculated asp þ ¼ ðp þ ρ m c mũ Þ=2 andp À ¼p Àp þ , respectively. This method is preferred to taking pressure time series at several locations, as in the experiment, as the CFD signal further downstream can potentially be expected to suffer from numerical dissipation to some extent. The up and downstream components are then processed in the same way as described for the experimental pressure measurements to calculate an impedance split into resistance and reactance.
The impedance can also be presented in non-dimensional form. Here the data is non-dimensionalised by using the magnitude of the impedance as ω-0. As there is no out-of-phase component of Z at ω ¼ 0 this will be the quasi-steady resistance. This can be found from the relation of mean pressure drop across the orifice to the ideal velocity (i.e. that predicted by the inviscid Bernoulli equation as found in the vena contracta) and hence to mean velocity U mean : where C d is the discharge coefficient, defined as The quasi-steady resistance is then Frequency is typically non-dimensionalised using a typical flow timescale to calculate the Strouhal number. Defined here as Here the mean velocity at the orifice, U mean ¼ C d ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2Δp=ρ 0 p , is obtained from CFD and experiment for the respective sets of results. By using the discharge coefficient from the steady-state CFD calculation errors in calculation of the mean flow field can be separated from those in the unsteady flow.

Mean flow field
The predicted mean axial velocity for the steady flow field is shown in Fig. 4. It can be seen that while the shortest hole shows a vena-contracta emerging from the orifice, the two longer holes show turbulent pipe flow behaviour in the holes before the exit into the plenum. This is confirmed by Fig. 5 which show the predicted streamlines entering the orifice. For L=D ¼ 0:5 the flow separates at the entrance and the flow does not reattach which gives the typical vena-contracta flow. Whereas the longer L=D ¼ 5 and 10 cases show that the flow reattaches after the initial separation.
The discharge coefficient, C d , can be found from both experiment and CFD using Eq. (25). The results are shown in Table 1. The largest discrepancy is between the predicted and measured values for L=D ¼ 0:5. Exact agreement for short holes is difficult to obtain as the C d will be determined by details such as the actual radius of curvature of the inlet to the orifice which cannot be matched in the CFD geometry. However the trend with L/D is well reproduced and the absolute agreement for the longer holes is good.

. Unsteady flow-field
The URANS formulation employed here allows for the unsteady flow field to be investigated computationally. Shown in Fig. 7 is the fluctuating axial velocity (u 0 ¼ u À u) at eight points in an acoustic cycle. The results are shown for the L=D ¼ 0:5 hole at a frequency of 420 Hz. To indicate the phase with respect to the pressure wave, the position in the cycle of the pressure at the entrance of the hole is indicated. As the signal amplitude is not high enough to reverse the flow in the orifice there are a succession of regions of increased and decreased velocity convected downstream from the hole. We can also see that the flow close to the wall in the orifice is the first to respond with the flow at the centre following later. A negative axial velocity fluctuation is seen at the edge of the orifice before the start of the negative Δp 0 phase of the cycle. This is most likely due to the increased recirculation (and hence negative axial momentum) seen in the detached region when the velocity in the centre is increased due to the positive Δp 0 part of the cycle. The process is reversed at the start of the positive Δp 0 part of the cycle when a positive axial velocity fluctuation is observed due to the decreased recirculation.
Reattachment (defined as the point where axial shear stress equals zero) occurs for the orifices of L=D ¼ 5 and 10. It is of interest to examine the movement of the reattachment point during one acoustic cycle. In Fig. 6 the variation of the length of the separation bubble is plotted against the cycle phase for 340 Hz. The phase angle is referenced to the pressure fluctuation at the upstream side of the orifice. It can be seen that the movement of the reattachment is almost in anti-phase to the pressure disturbance. The movement is small with respect to the mean distance of the reattachment point from the entrance of the orifice with this movement becoming smaller for the longer hole. It is thought that this motion will have negligible influence on the acoustically excited flow. However, the ability to investigate unsteady flow details of this sort, where experimental measurements cannot, is one of the advantages of this time-resolved CFD method. A quantitative comparison with the experimental impedance results is presented in the next section.

Orifice impedance data
The complex impedance of the short, L=D ¼ 0:5, orifice, expressed as resistance and reactance, is shown in Fig. 8. The impedance is shown as a function of frequency from both experiment and CFD calculation, with the impedance in both cases being calculated from the pressure measurement in the duct section. As described previously, at low frequencies the impedance is almost entirely resistive with the reactance being low. As frequency increases the reactance increases but does so at a faster than linear rate. For this orifice the resistance is seen to decrease as frequency increases. The experimental data is seen to be well reproduced by the CFD results.
For comparison we also include impedances calculated using the analytical models presented in Eqs. (9) and (13). The analytical models give reactances which are linear with frequency and are therefore not capable of reproducing the experimental trends seen here. More notably, the resistance predicted by the analytical models is almost constant with frequency and therefore give poor predictions at higher frequencies, highlighting the advantage of the CFD method. At very low frequencies (below 300 Hz) the Bellucci model is seen to give a closer prediction to the resistance than the CFD, however it should be noted that in this model the resistance was calculated using the experimentally measured discharge coefficient. This may not be available if the model is used as a design tool. Of particular interest is the ability of CFD to reproduce the negative resistance seen at high frequencies. This observation is supported by Fig. 9 which shows the reflection coefficient against Stouhal number (calculated as in Section 4.4) for this orifice. The reflection coefficient is defined as the ratio of the reflected wave amplitude in the duct to the incident wave amplitude which can be found as in Section 4.4. The experimental reflection coefficient goes above one at St ¼ 2:8 implying that the reflected wave is larger than the incident wave due to self excitation behaviour; i.e. some energy is being transferred from the mean to the acoustic part of the flow. Howe's energy corollary [39] shows that rotational kinetic energy can be converted into acoustic energy. For the shortest orifice we see in Fig. 5 that there is a significant region of recirculation in the orifice caused by separation. Fig. 7 shows that this recirculation is excited by the acoustic forcing which would lead to an increase in rotational kinetic energy which is then transferred back to the acoustic field. A similar

Strouhal Number
Experiment CFD Fig. 9. Reflection coefficient against Strouhal number from experiment and CFD for L=D ¼ 0:5. phenomenon has been observed in the work of da Silva et al. [16] where a recirculation at the open end of a pipe was seen to lead to energy being transferred to the acoustic field. It is particularly noteworthy that the URANS CFD method, capable of using with realistic geometries and Reynolds numbers, is able to reproduce this and further investigation of the CFD data to shed light on this self excitation is currently underway. The impedance for the longer holes is shown in Fig. 10 for L=D ¼ 5 and Fig. 11 for L=D ¼ 10. For both these cases the reactance again increases at a faster than linear rate with frequency which the analytical model of Eq. (13) is incapable of reproducing. The CFD method, however, captures the reactance very accurately. Unlike for the shorter orifice these cases show a very different behaviour for resistance as a function of frequency with the resistance increasing, rather than decreasing, with frequency as seen for the L=D ¼ 0:5 orifice. It is clear that a single analytical model of the sort seen in Eq. (13) will struggle to predict the very different trends for short and long holes whereas we see that the CFD approach is capable of correctly reproducing this trend. The absolute agreement between experiment and CFD for resistance is not as good at higher resistances. One possible explanation for this is that for the longer holes at higher frequencies the reflection coefficient increases. The impedance is given by wherep þ ð0Þ andp À ð0Þ are the up and downstream travelling pressure waves evaluated at the orifice plane, respectively. As the reflection coefficient increases the difference between the upstream and downstream travelling waves will reduce. Hence the impedance calculated as above will be susceptible to errors in measurement ofp þ andp À for both experiment and CFD. As the in-phase resistance becomes smaller relative to the out-of-phase reactance, it is likely the error will show up more for resistance than reactance. This may explain the non-monotonic nature of the experimental resistance results.

Frequency [Hz]
Resistance Experiment Resistance CFD Resistance Analytical Reactance Experiment Reactance CFD Reactance Analytical

Frequency [Hz]
Resistance Experiment Resistance CFD Resistance Analytical Reactance Experiment Reactance CFD Reactance Analytical Fig. 11. Impedance as a function of frequency for the L=D ¼ 10 orifice. Results from experiment and CFD. Also shown are impedances calculated using the analytical model of Eq. (13).
The ability of the CFD approach to correctly predict the trend of resistance and reactance with orifice length is shown in Figs. 12 and 13 which show resistance and reactance respectively for the three orifices. The resistance and reactance has been normalised as described in Section 4.4 and both are plotted against Strouhal number. When plotted together in this way the ability of the CFD to accurately predict the trend with orifice length can be seen. The agreement for reactance is good and for resistance the very different trend with Strouhal number for the short and long holes is captured. We also see that for the L/D ¼0.5 hole, normalising resistance and frequency as we have done here improves the agreement between experiment and CFD. The collapsing of the two sets of data onto one another here is mainly due to the effect of normalising the frequency to produce Strouhal number. This suggests that the negative resistance observed at higher frequencies for this hole is due to some self-excitation caused by an interaction of the acoustic and mean flow fields which occurs at a particular range of Strouhal number. For the longer holes, as discussed above, the resistance predictions may be affected by errors in measurement.

Conclusions
Experimental data for the specific impedance of plenum-fed orifices of L/D of 0.5, 5 and 10 has been obtained at a range of frequencies. A URANS CFD technique has been employed to reproduce this data. The CFD method used was an unsteady pressure-based compressible code with characteristic based boundary conditions used to give correct acoustic behaviour at the inlet and outlet of the computational domain. By making the plenum side of the domain sufficiently large, but still finite, and using a non-reflecting boundary it has been possible to reproduce the effect of the plenum which places a pressure node at the upstream entrance to the orifice. This plenum arrangement allows the acoustic response at a range of frequencies to be investigated using both experimental and CFD techniques.
The CFD prediction of reactance for all three holes was in good agreement with experimental data, both exhibiting a faster than linear increase in reactance with frequency which the analytical model cannot reproduce. The CFD method was found to correctly reproduce both the increase of resistance with frequency for L=D ¼ 5 and the decrease for L=D ¼ 0:5 which the analytical models were not able to do. The CFD was also able to capture the increase of reflection coefficient to values above unity at a Strouhal number of approximately 2.8 for the L=D ¼ 0:5 case. By normalising the resistance using the steady-state resistance based on C d , it was possible to separate the prediction of the unsteady behaviour from that of the mean flow-field, which for short orifices can be difficult to predict. By demonstrating this agreement with experiment it was shown that the URANS approach to turbulence modelling is capable of correctly predicting the acoustic response of a turbulent flow without the need to resolve in detail the process by which turbulence dissipates acoustically driven velocity fluctuations.
The numerical method demonstrated here offers two immediate avenues for further exploitation. While the CFD has been applied in this work to relatively simple geometries no tuning has been used to get the excellent agreement for the very different trends seen with the three geometries investigated. This, together with the use of an unstructured URANS CFD code, suggests that the method can be used to calculate the impedance of more complex aerodynamic flow features for which no analytical models are available, without resorting to the extra computational expense of LES or DNS methods. Secondly the three dimensional and time resolved nature of the results generated by this method mean that the details of such unsteady flows, for example the phase shift observed between the response at the centre and edge of the orifice, can be investigated in a way not possible in experiment.