On the scattering of entropy waves at sudden area expansions

In this work, we investigate both numerically and theoretically the sound generated by entropy waves passing through sudden area expansions. The numerical approach is based on a triple decomposition of the flow variables into a steady mean, a small-amplitude coherent part, and a stochastic turbulent part. The mean flow is obtained as the solution of the Reynolds-Averaged Navier-Stokes (RANS) equations. The equations governing the coherent perturbations are linearised and solved in the frequency domain. To account for the effect of turbulence on the coherent perturbations, a frozen eddy viscosity model is employed. When entropy fluctuations pass through the area expansion, the generated entropy noise behaves as a low-pass filter. The numerical predictions of the noise at low frequencies are compared to the predictions of compact, quasi-one-dimensional, and isentropic theory and large discrepancies are observed. An alternative model for the generated entropy noise tailored for area expansions is then proposed. Such model is based on the conservation of mass, momentum, and energy written in integral form. The model assumes zero frequency and the one-dimensionality of the flow variables far upstream and downstream of the expansion. The predictions of this model agree well with the numerical simulations across a range of finite subsonic Mach numbers including low, intermediate, and high Mach numbers.


Introduction
Small-amplitude perturbations superimposed on a single-component, compressible flow can be classified into three types of disturbances [1,2]: (i) acoustic waves that are homentropic and irrotational, (ii) entropy waves that are incompressible and irrotational, and (iii) vortical waves that are incompressible and homentropic. If the flow is uniform, and viscous effects and heat transfer are negligible, these waves are independent and can be considered separately. In this case, acoustic waves propagate in all directions in relation to the flow at the speed of sound, while entropy and vortical waves are simply swept by the flow silently.
Viscous terms, heat transfer, mean-flow gradients, and solid boundaries can act as coupling terms. For instance, acoustic waves can be scattered into vortical waves by sharp edges [3] or dissipated into heat (entropy) by viscous effects. Likewise, entropy waves can generate sound when accelerated/decelerated (mean-flow gradients) [4,5] or when interacting with solid boundaries, such as aerofoils [6][7][8]. This is termed entropy noise and is particularly relevant in the context of confined combustion, as occurs in the combustors of gas turbines [9] and rockets engines [10].
Unsteady combustion is, through unsteady heat release rate, a source of flow disturbances. In gas turbines, these disturbances are the origin of two distinctive sources of noise [11][12][13]: (i) direct and (ii) indirect sources. Unsteady heat release rate, through volumetric expansion [14], acts as a monopole source of acoustic fluctuations, the so-called direct combustion noise. Unsteady combustion also generates advective disturbances in the form of vorticity, temperature (entropy), and mixture-composition fluctuations. These are silent while advecting in the combustor but generate noise at the exhaust when interacting with the turbine blade rows [15]. This source of noise is termed indirect combustion noise and can be further decomposed as vortical [16][17][18], entropy [4,5] and, compositional noise [19][20][21] after the name of the flow perturbation which creates them. Entropy noise is believed to be its main component [9].
Entropy noise produced by the acceleration of entropy waves in nozzle flows was demonstrated experimentally by the Entropy Wave Generator (EWG) test rig at the German Aerospace Center (DLR) [5]. Leyko et al. [22] showed that the experimental results for sonic throat conditions can be reproduced using the compact theory of Marble and Candel [4]. This seminal work proposed an analytical solution of the linearised, quasi-one-dimensional, and isentropic Euler equations. The model assumes zero frequency and is valid for configurations where the wavelengths of both the entropy and acoustic waves are long compared to the length of the nozzle (the compact assumption).
Duran et al. [23] also compared the DLR EWG experimental results with the compact theory [4], in this case for subsonic conditions. They found that the model qualitatively predicts the shape of the experimental signals, but quantitative disagreement was observed. The mismatch was attributed to the noncompactness of the entropy waves for this condition owing to the lower convective velocity. Several efforts have been made to extend the frequency range of validity of the compact theory, including analytical solutions for nozzles with linear velocity profiles [24,25] and the use of asymptotic expansions of the governing equations in terms of frequency [26,27]. Duran and Moreau [28] proposed a solution based on the Magnus expansion method that is valid at any frequency within the range of validity of the quasi-one-dimensional assumption. Above a low-frequency limit, the entropy profiles become three-dimensional. This effect must be included in the equations of the acoustics that can be treated as onedimensional [29,30]. Further extensions of the quasi-one-dimensional theory include non-linear effects [31], compositional inhomogeneites [20,32], mean-flow heat trasfer [33], and viscous effects [34].
A limitation of the experimental campaign of the DLR EWG was the lack of measurements of entropy noise upstream of the nozzle (the so-called reflection coefficient). The study of this component is of great technical importance because it can modify the stability of the combustor [35,36] leading to largeamplitude, self-sustained thermoacoustic instabilities [37,38]. To address this limitation, a new rig was designed at the University of Cambridge, the Cambridge Wave Generator (CWG) [39], where the entropy spots were accelerated through an orifice plate. Rolland et al. [40] compared the experimental results with the predictions of quasi-one-dimensional and isentropic models [4,28], showing that the experimental reflection coefficient under subsonic conditions was significantly higher than predicted by such theories. The difference was attributed to the intrinsic differences between an isentropic nozzle, as described by the aforementioned theories, and an orifice plate, as used in the experiment. The orifice comprises a sudden area contraction followed by an expansion. In the contraction portion, the flow remains essentially attached and isentropic. However, in the expansion part, the flow separates, creating a large recirculation zone which leads to turbulent friction and strong mean-entropy gradients at the downstream side. De Domenico et al. [21,41] proposed the mean nonisentropicity produced at the area expansion as the source of the discrepancy and derived an alternative compact model that successfully reproduces the experimental results.
To obtain further physical insight into the entropy-noise generation mechanism in non-isentropic flows, Yang et al. [42] theoretically studied a sudden area expansion using an acoustic analogy. They found that, for low-Mach-number flows, the quasi-one-dimesional and isentropic theory of Marble and Candel [4] significantly overpredicts the actual reflection coefficient generated by entropy deceleration. The physical origin of the disagreement was explained in isentropic terms: the existence and length of the recirculation region creates a spatial gap between the source region and the geometrical area expansion and this produces a different acoustic response from when both are aligned. When the source term and the area expansion are aligned, the reflection coefficient was found to tend to the predictions of Marble and Candel [4] at low frequencies. The authors concluded that the main factor for the inability of the Marble and Candel model to capture the entropic response of sudden area expansions sustaining low-Machnumber bias flows is in fact the spatial delay of the source term associated with separation, and not the non-isentropicity of the flow as previously suggested.
The theoretical analysis of Yang et al. [42] was, however, restricted to low-Mach-number flows. In this work, we investigate both theoretically and numerically the sound generated by entropy fluctuations decelerated at sudden area expansions for several Mach numbers, including high subsonic Mach numbers.
The theoretical model proposed is an extension of the zero-frequency model developed by Ronneberger [43,44] for the acoustic scattering of area expansions. Here, we extend it to account for incoming entropy fluctuations. The approach only requires one-dimensionality of the flow variables at the inlet and at the outlet and places no dimensional restrictions in the close vicinity of the expansion itself.
The numerical approach is based on a linearisation of the Navier-Stokes equations about a turbulent mean flow. The linearised equations are then solved in the frequency domain to avoid spurious wave solutions that appear in the time domain [45]. A similar approach was used to successfully describe the acoustic response of perforations [46] and sudden area expansions [47,48] sustaining low-Mach-number bias flows. In these publications, two simplifications were considered: (i) the effect of turbulence on the acoustic field was neglected and (ii) the resolution of the energy equation was avoided by assuming that the mean flow was incompressible and the acoustics isentropic. The latter assumption is justified owe to the low Mach number of the configurations of interest. The first simplification was relaxed by [49,50] by decomposing the flow variables into mean, coherent and, turbulent parts (the so-called triple decomposition) as proposed by Reynolds and Hussain [51]. Additionally, some authors have considered the energy equation for laminar flows, e.g. Schulze et al. [52] to characterise perforations and Ullrich et al. [53] to investigate the acoustic-entropic coupling in a Laval nozzle. In this work, we consider both the effect of turbulence and the energy equation in the governing equations. To this end, we adapt the triple decomposition for compressible flows proposed by Alizard et al. [54] to the context of acoustic simulations.
The contributions of this paper are threefold. First, we present, for the first time in the context of acoustic simulations, a triple decomposition adapted to high-Mach-number, compressible flows. Secondly, we adapt a compact theory tailored for area expansions to account for incoming entropy fluctuations. Thirdly, we use these two approaches to systematically characterise the unsteady response of a subsonic area expansion forced by incoming entropy fluctuations, a configuration that represents a canonical example of internal flows with largescale separation and stagnation pressure losses (non-isentropicity).
This article is structured as follows. First, the configuration of interest is introduced in Section 2. Secondly, the numerical and analytical approaches are presented in Sections 3 and 4, respectively. Results are then given in Sections 5 and 6 for the mean flow and the perturbation response, respectively. Finally, conclusions are drawn in Section 7.

Problem formulation
We consider a canonical area expansion consisting of two concentric cylinders with different cross-sectional areas as sketched in figure 1. A subsonic, uniform flow in the streamwise direction is imposed on the left-hand side of the domain. The left-hand-and right-hand-side ducts are denoted, respectively, as upstream and downstream ducts and all the variables defined in them are denoted by the subscripts (·) u and (·) d , respectively. The expansion ratio, η = A u /A d , defined as the ratio of the cross-sectional areas of the upstream, A u , and downstream ducts, A d , fully determines the geometry of the problem. The backplate, defined as the portion of the external wall located where the upstream duct expands, has an area The variables defined on the backplate are denoted by the subscript (·) b . The geometry exhibits cylindrical symmetry around the central axis and, thus, can be described using a cylindrical system of coordinates centred in the axis at the location of the expansion. The axial, radial, and azimuthal coordinates are denoted by x, r, and θ, respectively, and its unit vectors by e x , e r , and e θ , respectively. Additionally, we assume that the flow variables are axisymmetric.
We assume that the flow at the inlet consists of a uniform component of velocityū u =ū u e x , densityρ u , and speed of soundc u , on which there is superimposed a small unsteady motion. Consequently, the flow variables can be decomposed into a steady mean, denoted by(·), and a perturbation component, denoted by(·). The mean flow velocity is normalised by the speed of sound to define the Mach number, i.e. M = |ū|/c. The perturbation component is small and, thus, the equations governing its evolution can be linearised around the mean flow allowing us to use the harmonic ansatz where ω is the angular frequency, t the time, and i the imaginary unit. The angular frequency can be normalised by the diameter of the upstream duct, D u , and by either the mean velocity in the upstream duct,ū u , to define a Strouhal number or by the speed of sound in the upstream duct,c u , to define a Helmholtz number Far from the area expansion, the acoustic field can be assumed to be composed of plane waves propagating in the same and opposite direction of the mean flow and denoted by(·) + and(·) − , respectively. The perturbation pressure and axial velocity can be written in terms of the plane waves asp =p + +p − and u = (p + −p − ) / (ρc) , respectively.
The plane acoustic waves can be classified into incoming waves, which are inputs to the system, and outgoing waves, which are outputs. Here, we consider uniform entropy waves imposed at the inlet,ŝ + u , as the input, with the two outputs being the reflected,p − u , and transmitted,p + d , acoustic waves. The incoming entropy fluctuations are normalised by the specific heat at constant pressure c p , such that σ =ŝ + u /c p . The two outgoing acoustic waves are normalised as where γ denotes the adiabatic index. The reflection and transmission coefficients are then defined as P − u /σ and P + d /σ, respectively. In this work, we seek to retrieve these reflection and transmission coefficients using two different techniques: a numerical approach, presented in Sec. 3, and an analytical model, presented in Sec. 4.

Numerical modelling
In this section, we present the numerical approach used to characterise the acoustic response of the sudden area expansion. The derivation of the governing equations starts from the compressible Navier-Stokes equations. A triple decomposition is then introduced (Sec. 3.1) and equations for the mean and perturbation parts are derived in Sec. 3.2 and Sec. 3.3, respectively.

Governing equations
The working fluid is taken to be a compressible, calorically perfect gas. The equations describing the flow are the conservation of mass, momentum, and energy that read, respectively, ∂ρ ∂t + ∇ · (ρu) = 0, (5a) where ρ is the density, u the velocity, p the pressure, and e t the total energy. Heat radiation is neglected. The fluid is taken to be Newtonian and, thus, the viscous stress tensor, τ, is given by where µ is the molecular dynamic viscosity and δ is the Kronecker delta unit tensor. Using the Fourier law, the heat-conduction vector, q, is defined as where T is the temperature and P r the Prandtl number. The temperature is related to the total energy by the relation e t = c v T + u · u/2, where c v is the specific heat at constant volume. Additionally, the equation of state relates temperature, density and pressure as p = ρRT, with R the specific gas constant. As explained in Sec. 2, the flow variables are decomposed into mean and perturbation parts. Following Reynolds and Hussain [51], such decomposition can be extended to account for turbulent motions using the so-called triple decomposition. The triple decomposition of the density and pressure fields into mean flow, coherent perturbations, and small-scale random motions reads with the definition of each component given bȳ The phase-average operator, · , is defined by with T the period of the harmonic motion, i.e T = 2π/ω. Note that the coherent component contains the acoustic motions, as well as vortical and entropic fluctuations. When considering compressible flows, it is convenient to introduce a densityweighted triple decomposition as proposed by Alizard et al. [54]. This decomposition is applied to the velocity, u, total energy, e t , and temperature, T , fields, such that where the subindex (·) f denotes density-weighted components (also called Favreaveraged components). The definitions of these terms arē To simplify the notation and to keep the parallelism with the laminar formulation of the linearised Navier-Stokes equations, the subindex (·) f is dropped hereafter. These decompositions are now introduced in the conservation equations, Eq. (5), to obtain the governing equations for the mean flow (Sec. 3.2) and coherent perturbations (Sec. 3.3).

Mean flow
Taking the time average of Eq. (5) and neglecting terms higher than secondorder of the coherent component, we obtain the governing equations for the mean flow variables, such that τ andq are the Favre-averaged stress tensor and heat conduction vector, respectively. Note that Eq. (13) is not exact and additional assumptions are required for its derivation (see [55] for further details). Eq. (14) contains two terms that need closure models: the Reynolds stress tensor, −ρu ′ ⊗ u ′ , and the turbulent heat flux vector, −ρu ′ c p T ′ . The Boussinesq approximation [56] is used to model the Reynolds stress tensor by assuming that it is proportional to the trace-less mean strain-rate tensor, i. e.
where µ t is a scalar termed eddy viscosity and k is the turbulent kinetic energy, i.e. k = 1 2 ρu ′ · u ′ . The turbulent heat flux vector is modelled using a Reynolds analogy, such where P r t is the turbulent Prandtl number. The governing equations are solved using the open-source finite volume solver OpenFOAM (version 4.1) [57]. A steady state solution is obtained using the SIMPLE algorithm [58]. The eddy viscosity, µ t , is determined using the komega-SST model [59].

Coherent perturbations
Taking the phase average of Eq. (5), subtracting Eq. (13), and neglecting terms higher than second-order of the coherent component, we obtain the governing equations for the coherent perturbations as follows where : is the double dot product, such that M : N = tr(M N ⊤ ). The fluctuation total stress tensor is given bỹ The fluctuating total heat flux vector is given bỹ The perturbation temperature is related to pressure and density through the linearised equation of state, that reads These equations (Eq. (17)) are termed hereafter the coupled formulation and are exact (within the assumptions given so far). To investigate the effect of mean non-isentropicity on the perturbation response of the area expansion, this exact formulation (Eq. (17)) is further simplified in the next section.

Uncoupled formulation
In this section, we simplify the perturbation equations by neglecting meanentropy gradients. To this end, the energy equation (Eq. (17c)) is given in terms of entropy as For common acoustic applications [46][47][48], the generation of entropy due to viscous effects and thermal conduction can be neglected in the derivation of the perturbation equations. Additionally, we can assume that entropy fluctuations due to turbulence are negligible, i.e. s =s f +s f . Taking now the phase average · of Eq. (21) and subtracting its mean, we obtain, after linearisation, the following equationDs where the operatorD/Dt = ∂/∂t +ū · ∇ denotes the convective derivative associated with the mean flow. Again, the subindex (·) f (denoting densityweighted variables) has been dropped to simplify notation. The first term of this equation shows that entropy is convected by the mean flow. The second term represents a source of entropy due to the interaction of acoustic and/or vorticity waves with mean-flow entropy gradients. For certain conditions (including low-Mach-number flows [46][47][48]), Eq. (22) can be further simplified by neglecting the mean-entropy gradient, i.e. ∇s ≈ 0, leading to a convective equation, such that This result effectively decouples the energy equation from the mass and momentum conservation equations. Eq. (23) together with Eqs. (17a) and (17b) is thus termed the uncoupled formulation. The relation between density, pressure, and entropy can now be obtained by linearising the Gibbs relation. The resulting expression is then integrated from a reference state and, subsequently, the phase average of that equation is taken. Neglecting also temperature fluctuations due to turbulence, i.e. T ′ f = 0, we obtainp This expression is the generalisation of previous work [47,49,50] to include entropy fluctuations.
For low-Mach-number flows, variations of the mean densityρ are negligible. Assuming also that turbulence fluctuations do not generate significant density variations, i.e. ρ ′ = 0, as in previous studies [49,50], we can prove that the density-weighted variables given by Eq. (12) simplify to the variables defined in Eq. (9) for the standard triple decomposition. This proves that this formalism collapses to the formulation proposed previously in the literature [49,50] for low-Mach-number flows.

Closure model
Following [51,54], and analogous to the closure model used for the mean flow, the fluctuations of the Reynolds stress tensor are assumed proportional to the fluctuations of the strain-rate tensor as Additionally, we assume that the turbulence is not affected by the acoustic fluctuations which allows the eddy viscosity to be considered "frozen" and equal to the value obtained for the mean flow. The turbulent heat flux vector (or dissipation of internal energy) is modelled again using a Reynolds analogy, so that

Finite element strategy
Eqs. (17a), (17b), and (17c), are referred as the coupled formulation while Eqs. (17a), (17b), and (23) are termed the uncoupled formulation. The two sets of equations are recast in the frequency domain and solved using a finite element method stabilised with a least-squares formulation [60]. The equations are implemented in the open-source computing framework FEniCS [61,62]. After discretisation, the resulting linear system is inverted using the sparse direct solver MUMPS [63,64]. The numerical implementation has been validated with several canonical test cases, including an isentropic nozzle (see Appendix A) and the acoustic scattering of an area expansion (Appendix B). Figure 2: Schematic of the region V used in the derivation of Eq. (27). S is the surface bounding V: At the inlet and outlet of the domain, non-reflecting boundary conditions are imposed. The implementation of the boundary conditions is based on the splitting of the flux vector defined at the boundaries (that appears naturally in the formulation after integration by parts) into incoming and outgoing waves in a similar way as that of the characteristic method [65]. This method allows incoming acoustic/entropy waves to be easily imposed at both boundaries while assuring non-reflectivity of the outgoing waves.

Analytical modelling
Let us define the control volume V bounded by the surface S as depicted in Fig. 2. The sub-surfaces S u and S d are taken sufficiently far from the expansion so that the flow variables can be considered one-dimensional over them. The conservation equations, given by Eq. (5), are integrated over this volume as the starting point for the derivation of the model. The volume terms are then transformed to surface integrals using the divergence theorem. The flow variables are assumed to evolve quasi-steadily, so that ∂(·)/∂t = 0. This assumption restricts the model to low frequencies. Viscous friction on the wall boundaries (S wall ∪ S b ) is neglected. Viscous effects, however, are implicitly retained in the rest of the domain. Finally, a non-slip and adiabatic boundary condition is imposed on the wall boundaries. Considering symmetry arguments, the conservation of mass, momentum, and energy, can be written [44], respectively, as where h t denotes total enthalpy. The relation between total energy and enthalpy, h t = e t + p/ρ, was used to obtain Eq. (27c). The pressure on the backplate, p b , in the momentum equation (Eq. (27b)), requires closure. A common assumption for separated flows [44,[66][67][68] is that the pressure acting on the backplate is equal to the pressure just before the area expansion, i.e. p b = p u .
The decomposition introduced in Sec. 2 is now used in Eq. (27) to obtain the governing equations for the mean (Sec. 4.1) and perturbation (Sec. 4.2) components.

Mean flow
The governing equations for the mean flow [69] are where β and Ξ are non-dimensional parameters defined as β =ρ d /ρ u and Ξ =c d /c u . Eqs. (28) form a system of three independent, non-linear equations with six unknowns. This system can be solved numerically after imposing three of them. Here, we impose the geometry through η, the type of gas through γ, and the normalised mass flow rate through the Mach number at the inlet, M u . The mean pressure ratio, defined as α =p d /p u = βΞ 2 , could be used alternatively in Eq. (28). The degree of non-isentropicity of the area expansion is described by the mean-entropy jump across the area expansion, given by

Perturbations
The governing equations for the perturbation part are now obtained by linearising Eq. (27). The resulting system of equations is rewritten in terms of plane waves propagating in the same and opposite direction to the mean flow,p + ,p − andŝ + , and arranged into input and output vectors to the system, written, respectively, as The resulting system of equations is given by and For low Mach number flows, i.e. M 2 u , M 2 d ≪ 1, the solution to Eq. (28) is Ξ = 1, β = 1 and M d = M u η. At this limit, the acoustic response when an entropy wave enters the system is identically zero, i.e.p − u = 0 andp + d = 0.

Mean-flow results
In this section, we report the mean-flow results obtained for a canonical area expansion as introduced in Sec. 2 using both the numerical and analytical approaches outlined in Sec. 3.2 and Sec. 4.1, respectively.
The mean-flow variables are determined by five non-dimensional parameters: the expansion ratio η, the Reynolds number Re = (ρ uūu D u ) /µ, the adiabatic index γ, the Mach number at the inlet M u , and the Prandtl number P r.
We consider the working fluid to be air, which fixes two of the parameters: γ = 1.4 and P r = 0.7. The expansion ratio selected in this paper is η = 0.346, due to the wealth of numerical [47,48,70] and experimental [44] data available for this geometry. The Reynolds number is taken to be Re = 10 6 , which is sufficiently high to assure that the flow topology is independent to this parameter (see [71]). Finally, three upstream Mach numbers are investigated throughout this publication: M u = 0.2, 0.5, and 0.8. The range of Mach numbers considered in this study includes subsonic high Mach numbers, for which few previous studies exist.
The upstream duct is modelled as a laminar, uniform flow. To achieve this, we impose uniform velocity and temperature profiles at the inlet, so that the targeted Mach number, M u , is obtained. The turbulent kinetic energy at the inlet is set to zero, i.e. k = 0. A zero-gradient condition is prescribed to both the pressure and specific rate of turbulent dissipation (omega). A slip boundary condition, that readsū ·n = 0, is imposed on the wall of the upstream duct. A non-slip, adiabatic boundary condition is imposed on the walls at the backplate and downstream duct. The viscous sublayer is fully resolved on the wall boundaries, which avoids the use of wall functions. At the outlet of the domain, uniform pressure is imposed and zero-gradient for the rest of variables.
Due to the axisymmetry of both the domain and boundary conditions, the mean flow is assumed axisymmetric. To exploit this, the flow is solved in a wedge domain of angle 5 • with appropriate boundary conditions in the azimuthal direction. The mesh extends only one cell in the azimuthal direction. The upstream and downstream ducts extends 50D u each in the streamwise direction. This length was selected to assure that an accurate separation of acoustic waves into downstream/upstream-propagating waves was possible at low frequencies. All the meshes are structured and composed mostly of hexahedral cells together with a layer of prismatic cells in the axis of revolution. A total of 276, 400 cells was used for each configuration. The meshes used yielded y + values around 1 at the duct walls. In the central part of the domain, a jet is formed. This jet remains almost parallel just after the end of the upstream duct. This translates to the pressure just before the expansion being equal to the pressure acting on the backplate as shown in figure 3(c). The central jet is separated from the low-speed recirculation zone by an expanding shear layer. At the end of the recirculation region, this shear layer merges with the boundary layer on the duct wall. After that point, the velocities at the centre and periphery of the duct slowly converge forming a fully-developed turbulent pipe flow at about 20D u downstream of the  expansion. This whole process creates strong turbulence as visualised in the form of eddy viscosity in figure 3(d).
The strong turbulent motion across the shear layer dissipates kinetic energy into heat by friction effects. This leads to a sharp increase of temperature and entropy across it as depicted in figure 4. The variation of entropy, which is almost constant in the upstream and downstream ducts, occurs in a very thin region and can be thought as a jump. This region will act as a source of entropy perturbations in the presence of acoustic/vortical waves, as described by Eq. (22).
The performance of the model is now assessed against the numerical results. To this end, the flow variables in the upstream and downstream ducts are sampled at the cross-sectional locations x = −1 and x = 30, respectively, and averaged over the surfaces. In the upstream duct, the flow variables are visually one-dimensional right up until the expansion and variations in the axial  direction are negligible as observed in figures 3 and 4. This justifies the selection of the sampling location very close to the area expansion. Figure 5 shows that the agreement between the numerical results and the model predictions is excellent for all the range of subsonic Mach numbers. Figure 5 also shows that, at low Mach numbers, the strength of the entropy jump is low and, thus, the values of the mean-flow variables are close to the predictions for an insentropic expansion. In contrast, when the Mach number increases, the strength of the entropy jump rapidly increases and the results deviate from isentropic theory. The increase of entropy is in fact directly linked to losses of stagnation pressure. This is observed in figure 5(c) that shows that, at high Mach numbers, the pressure recovered along the recirculation region is significantly lower for the sudden area expansion than for an isentropic expansion.

Scattering of entropy waves
In this section, the acoustic response of the expansion to incoming entropy waves is investigated using the numerical and analytical approaches introduced in Sec. 3.3 and Sec. 4.2, respectively.
We first characterise the expansion numerically, for which the mean-flow variables presented in Sec. 5 are interpolated to an acoustic mesh. This mesh is two-dimensional, fully unstructured and composed of 558, 339 triangles. The upstream and downstream ducts are 15D u and 40D u long, respectively. The coupled and decoupled formulations (see Sec. 3.3.3) of the linearised equations are formulated in cylindrical coordinates and solved assuming axisymmetry of the flow variables. The numerical implementation allows for the order of the basis functions to be selected from first-to fifth-order Lagrangian polynomials. The results presented here were obtained using second-order polynomials.
An incoming entropy perturbation of normalised amplitude σ =ŝ + u /c p is imposed at the inlet of the domain. Non-reflecting boundary conditions are prescribed at both ends, so only reflected and transmitted plane acoustic waves exist at the far upstream and downstream ducts, respectively. At the backplate, a rigid wall non-slip boundary condition is imposed, namelyũ = 0. Following [47], a slip boundary condition, that readsũ · n = 0, is prescribed on the walls of the upstream and downstream ducts. This approach reduces the computational cost as the acoustic boundary layer does not need to be resolved. However, any viscous dissipation associated with the walls is not captured. Gikadi et al. [48] showed that, for this geometry, the acoustic results were the same using either slip or non-slip boundary conditions on the duct walls.
We consider frequencies ranging 8 × 10 −4 ≤ He ≤ 1.2 for M u = 0.2, and 8 × 10 −4 ≤ He ≤ 3.0 for M u = 0.5 and M u = 0.8. In the upstream duct, the shortest acoustic wavelength is λ ≈ 0.42 obtained for M u = 0.8 and He = 3.0 and for which the resolution is approximately 10 elements per wavelength. In the downstream duct, the shortest acoustic wavelength is λ ≈ 2.5 obtained for M u = 0.5 and He = 3.0, which yields a resolution of around 63 elements per wavelength at the coarsest section of the duct. The shortest wavelength for the entropy at the upstream duct is λ ≈ 1.04 obtained for He = 3.0 and M u = 0.5. There are approximately 21 elements per entropy wavelength at this frequency. A detailed analysis of the error obtained for a given number of points per wavelength and polynomial degree for the current implementation can be found in Appendix C.
Figures 6(a) and 6(b) show the spatial distribution of entropy obtained using the coupled formulation of the linearised Navier-Stokes equations (LNSE). Neglecting viscous dissipation and heat conduction, the evolution of the perturbation entropy is given by Eq. (22). This equation comprises two terms: a convective term, that describes the advection of entropy by the mean flow as a passive scalar, and a source term that describes the generation of perturbation entropy through the interaction of acoustic/vorticity waves with mean-entropy gradients. At low Mach numbers, figure 5(b) shows that the mean-entropy jump for this configuration is weak and, thus, the source term in the aforementioned equation is negligible. In this case, the perturbation entropy is simply advected by the mean flow as observed in figures 6(a) and 6(b) for two different frequencies. For both frequencies, the entropy presence expands after the expansion, but does so gradually, meaning that there is an annular region in the recirculation zone where the entropy is negligible. At higher Mach numbers, the source term becomes significant and further perturbation entropy is generated after the expansion in the thin region where the mean entropy sharply increases (see figure 4(b)).
Yang et al. [42] showed that, at low frequencies and low Mach numbers, the main acoustic source term due to entropy is a dipole given by The mean pressure,p is constant upstream of the area expansion, it then recovers along the recirculation region -in the streamwise direction -and becomes constant again after it, as shown in figure 3 (c). The entropy perturbation,s, and, hence, the entropy-related source term (Eq. (31)) are only significant within the expanding jet, with negligible presence in the recirculation zone. Figure 6(c) shows an example of acoustic field generated by the entropy fluctuations.
To compute the reflection, P − u /σ, and transmission, P + d /σ, coefficients, the magnitude of the acoustic plane waves is extracted using the multi-microphone method [72,73]. To this end, two post-processing zones are selected: −15 ≤ x ≤ −5 upstream and 25 ≤ x ≤ 38 downstream of the area expansion. The variables are averaged over the cross-sectional area at each location to reduce the influence of numerical disturbances. 201 and 131 equispaced points are used in the upstream and downstream ducts, respectively. Figures 7, 8, and 9 show the magnitude of the reflection and transmission coefficients for M u = 0.2, M u = 0.5, and M u = 0.8, respectively. The reflection coefficient behaves as a low-pass filter: at low frequencies, the reflection amplitude is nearly constant and, when the frequency increases non-compact effects cause a drop of its value. The drop-off can be explained by Eq. (31) and figures 6(a) and 6(b): when more than a half period of the entropy wave is contained within the source region, cancelling effects between positive and negative parts of the wave occur. As for the transmission coefficient, its amplitude remains nearly constant at low frequencies, then a bump appears at moderate frequencies to finally drop at higher frequencies. The ratio between the maximum value at the peak and the value at low frequencies increases with increasing Mach numbers. For high Mach numbers, the value of the peak almost doubles the value at compact frequencies.
At low frequencies, the agreement of the model predictions with the numerical results (coupled formulation) is excellent for the three Mach numbers. Note that the coupled equations are an exact formulation (within the general assumptions given in Sec. 3.3) of the problem and the results obtained with it should be used as reference. The uncoupled formulation, on the other hand, neglects the mean-entropy gradient and its predictions should be used to quantify the validity of that assumption. For M u = 0.2, the results obtained for the coupled and uncoupled solvers are very close. Some discrepancies between the formulations are observed for the transmission coefficient, that can be attributed to postprocessing errors in the computation of the magnitude of the coefficients. This error is produced by a combination of: (i) the very long acoustic wavelengths that turns the application of the multi-microphone method into an ill-posed problem and (ii) the very low values of the transmission coefficient and, consequently, the low ratio between the actual signal and noise, either numerical (introduced by instance by the interpolation of the mean flow) or vortical waves. When the Mach number increases, the results of the coupled and decoupled solvers show some discrepancies. The reason is that the decoupled formulation neglects the perturbation entropy source related to the mean-entropy gradient and, as observed in figure 5(b), this becomes significant when increasing the Mach number.
In figure 10, the numerical results obtained at low frequencies are compared to the predictions of the compact model derived by Marble and Candel [4] for ducts with smooth area variations, i.e. nozzle flows. Because this model is compact, its predictions only depend on the conditions upstream and downstream of the area expansion. If the assumptions used for its derivation are valid, it should be a suitable candidate for sudden area expansions as considered here. These  assumptions are that the flow is inviscid, quasi-one-dimensional, and isentropic. For its derivation, Marbel and Candel imposed the conservation of mass, stagnation enthalpy, and entropy. The value of the reflection coefficient predicted by that model is higher than the actual value across the whole range of Mach numbers. The mismatch is significant at low Mach numbers (for M u = 0.2 the relative error is approximately 117%) and it becomes very severe when the Mach number is large: Marble and Candel predictions are two orders of magnitude higher that the predictions of the proposed model when the Mach number is close to unity. The predictions of the transmission coefficient also differ significantly, but the values remain of the same order of magnitude. The mismatch slightly increases with increasing Mach number: the relative error for M u = 0.2 is approximately 94% and increases to 121% for M u = 0.8.
The strong discrepancy between the predictions of the quasi-1D, isentropic model of Marble and Candel and the actual computational results was explained by Yang et al. [42] for low-Mach-number area expansions. They showed that the entropy-related acoustic response has a very strong dependence on the length of the recirculation region: when the flow reattaches immediately after the separation, i.e. the geometrical expansion matches the flow expansion, the values of the acoustic coefficients tend to the predictions of the Marble and Candel model. However, they significantly deviate when the length of the recirculation zone is increased. There exists a length after which an increase of such length does not bring any changes to the values of the acoustic coefficients. In terms of acoustic sources, the previous observation states that the Marble and Candel model is valid when the source is located at the actual expansion, but not when the entropy-related source is moved to a different position, in this case, into the downstream duct.
The misalignment of the geometry and the acoustic source is sufficient to explain the discrepancies observed at low Mach numbers. However, when the Mach number increases, the mean-entropy gradient increases, and, thus, there exists a source of fluctuation entropy within the domain (see Eq. (22)). In other words, the isentropicity assumption used in the Marble and Candel model is no longer valid. To explore the importance of this assumption, we slightly modify the model proposed in Sec. 4. The model relies on the mean flow variables obtained using the equations presented in Sec. 4.1. The perturbation variables are obtained by imposing the linearised conservation of mass and momentum (Eqs. (27a) and (27b)), but instead of the conservation of stagnation enthalpy, we impose the conservation of perturbation entropy, that readss u =s d . This model captures the source/geometry mismatch through the momentum equation and the hypothesisp u =p d , but also imposes the linear isentropicity condition. Figure 11 shows the predictions of the model. We observe that the model predicts values close to the results obtained at low frequencies with the uncoupled numerical solver. This is expected, since both approaches neglect perturbation entropy sources. The values of both the reflection and transmission coefficients at low and intermediate Mach numbers obtained when neglecting the sources of perturbation entropy are close to the actual values (the results of the coupled solver). When the Mach number goes beyond M u ∼ 0.5, the values diverge. The disagreement is especially large for the transmission coefficient: the actual value monotonically increases until unitary Mach, while the isentropic model drops after M u ≈ 0.7 reaching a zero value when the Mach number reaches one. As for the reflection coefficient, the constant-entropy model slightly overpredicts the actual value (a relative error of about 18% is obtained for M u = 0.8), but the general trend is correctly predicted. These results show that at low and moderate Mach numbers, the geometry/source misalignment is the dominant factor to explain the strong divergence between the Marble and Candel model and the numerical predictions obtained for area expansions. At high Mach numbers, non-isentropic effects become significant: for the reflection coefficient, perturbation entropy sources are still a second order effect, but they highly impact the transmission coefficient.
Finally, we present a parametric study performed with the model presented in Sec. 4. Figure 12

Conclusions
The acoustic response of a generic area expansion subjected to incoming entropy waves was investigated both numerically and analytically. The numerical approach is based on a triple decomposition of the flow variables into a steady mean, a coherent perturbation part (containing acoustic, vortical, and entropy waves), and a stochastic turbulent part. Separated governing equations for each part were derived. The mean flow was obtained as the solution of the Reynolds-Averaged Navier-Stokes (RANS) equations closed by the k-omega-SST turbulence model. The amplitude of the coherent perturbations was assumed to be small and the equations governing this part were linearised and solved in the frequency domain. To account for the effect of turbulence on the coherent perturbation, a frozen eddy viscosity model was used. From a post-processing point of view, the use of the eddy viscosity eases the computation of the acoustic waves because it damps the vortical waves generated at the expansion. Two different formulations of the perturbation equations were adopted in this work: (i) an exact formulation, that was termed a coupled solver, and (ii) a formulation that decouples the energy equation from the mass and momentum conservation equations by neglecting mean-entropy gradients, termed an uncoupled solver.
The results show that the area expansion behaves as a low pass-filter to entropy fluctuations. The reflection coefficient has a virtually constant value at low frequencies and drops at higher frequencies due to non-compact effects. The transmission coefficient, on the other hand, also has a constant value at low frequencies, but increases its value at intermediate frequencies before dropping at higher frequencies. This behaviour was found to be consistent across the range of subsonic Mach numbers studied here. The values of the coefficients at low frequencies were found to be significantly over-predicted by the quasi-onedimensional and isentropic theory of Marble and Candel [4]. The mismatch is especially pronounced for the reflection coefficient and for subsonic high Mach numbers (for M u = 0.8 the predictions of the model are one order of magnitude higher than the actual values).
Two main physical mechanisms can be responsible for the disagreement: (i) the mean-flow non-isentropicity (in other words, stagnation pressure losses) [41] or the spatial misalignment between the source term and the geometrical area expansion owe to the presence of the recirculation zone [42]. To explore the relative importance of each of them, numerical simulations were performed with the decoupled solver. This solver neglects the mean-entropy gradient term in the energy and, hence, any source of perturbation entropy. For the reflection coefficient, it was found that, for low to moderate Mach numbers, the predictions of the coupled and uncoupled solvers were similar. Some disagreement was observed for higher Mach numbers (M u 0.7), but not enough to explain the discrepancies with Marble and Candel. These results suggest that this model fails to predict the reflection coefficient of area expansions mainly owe to the spatial delay of the source term associated with separation. For the transmission coefficient, on the other hand, the trends predicted for both solvers are similar up to M u ≈ 0.6, where they diverge. After that point, the errors of the Marble and Candel model and the uncoupled solver are similar. This suggests, that for low to moderate Mach numbers, the main factor of discrepancy is again the spatial delay associated with separation, but beyond that point the mean flow non-isentropicity is as important.
An alternative compact model was developed based on the conservation of mass, momentum, and stagnation enthalpy. The model captures the threedimensional effects through the momentum equation. The agreement between the numerical results and the model predictions is excellent for all the range of subsonic Mach numbers studied. Finally, a parametric study of the acousticentropic response of the area expansion was performed based on the model. It was found that the reflection coefficient for large area expansions (η 0.2) tends to zero for any subsonic Mach number. This is in clear opposition to divergent nozzles with large area expansions, where this coefficient monotonically increases with decreasing expansion ratios. This result has important implications when modelling entropy acceleration through orifices using two-step models for perforations as showed by De Domenico et al. [21,41] Forcing The mean flow is computed using a steady-state solver for the non-linear compressible Euler equations in conservation variables. The equations are discretised in space using a continuous-Galerkin formulation stabilised using the least-squares method [60]. The discretised non-linear problem is solved using a fully-implicit, pseudo-time-stepping algorithm [74]. Second-order interpolation polynomials (p = 2) were used for the computations. The flow, which is fully defined by a Mach number at the inlet of M u = 0.0212, is accelerated at the nozzle to a Mach number of M d = 0.7029 at the downstream duct, as observed in figure A.13.
To compute the perturbations variables, both the molecular and turbulent viscosities are set to zero in Eq. (17) and the coupled equations are solved in cylindrical coordinates assuming axisymmetry of the flow variables. The equations are solved for a normalised frequency of He = 0.004 and for three types of incoming excitations: (a) downstream-propagating acoustic waves at the inlet P + u , (b) upstream-propagating acoustic waves at the outlet P − d , and (c) entropy waves at the inlet σ + u . Again, second-order interpolation polynomials (p = 2) were used for the computations. The results are decomposed into upstream/downstream propagating components using the multimicrophone method in two post-processing zones spanning −200 ≤ x/R d ≤ −50 and 50 ≤ x/R d ≤ 200 for the upstream and downstream ducts, respectively. The results obtained are summarised in table A.1, together with the predictions of the Marble and Candel model. The agreement is excellent for the three types of forcing, with the values being in agreement to, at least, the third decimal place.

Appendix B. Validation case: Acoustic response of the area expansion
To further validate the code, we compare the numerical results obtained with the present implementation with the experimental dataset for an area expansion of Ronneberger [44]. This experiment has been used to benchmark several numerical approaches [46,48,70,75]. Most of the parameters used in the experiment are the same as the ones used throughout this paper: the expansion ratio, η, the working fluid and the Mach number at the inlet M u = 0.2. The Reynolds number of the experiment corresponds to Re ≈ 10 5 which is lower than the value used in this publication Re = 10 6 . A fundamental difference of the results presented in this paper with the experimental measurements is that the latter has a boundary layer at the wall of the upstream duct, while here we assume a uniform profile. In order to perform a fair comparison, an additional mean flow is computed prescribing a non-slip boundary condition at the wall of the upstream duct. At the inlet, a fully developed velocity profile is prescribed. Such profile was obtained from an additional incompressible simulation of a uniform, straight duct. The turbulent variables of the developed profile are also specified. Additionally, a zero-gradient is imposed to the pressure and a uniform profile to the temperature. The domain is the same described in Sec. 3.2. A fully structured mesh composed of 795, 000 cells is used here. Again, the viscous sublayer was resolved and, thus, the mesh was designed so that y + ∼ 1 is satisfied at the walls. Two Mach numbers were simulated: M u = 0.2 and M d = 0.5. Figure B.14 compares the results with the uniform case showing good agreement.
The acoustic response is then computed using both mean flows (slip and non-slip upstream wall) for M u = 0.2. The boundary conditions, meshes and numerical parameters are defined in Sec. 6. The simulations were performed using the coupled solver. The relations between acoustic waves are summarized in the scattering matrix [76] which is defined as (B.1) Figure B.15 shows the results obtained for the scattering matrix. As can be seen, the simulated results agree well with experimental data using both mean flow fields. T − is underpredicted, but this trend is consistent with previous numerical studies [47,70,75]. For higher frequencies, the acoustic results obtained for separations with a uniform flow and a flow with a boundary layer in the upstream duct diverge.

Appendix C. Numerical accuracy
In this appendix, we study the accuracy of the numerical implementation proposed in this paper for a given number of points per wavelength and order of the interpolation polynomials. To this end, we solve the perturbation equations for the coupled formalism (Eq. (17)) in a fixed domain [0, 1] × [0, 1]. The domain is discretised by five fully non-structured meshes composed of triangles. Each mesh has the same number of elements prescribed to each of its sides and this number is used to define the number of elements per wavelength used in figure C. 16.
Non-reflecting boundary conditions for the acoustic waves are prescribed at left-hand and right-hand sides of the domain. Additionally, a left-to-rightpropagating acoustic wave is introduced at the left side. The variables are assumed axisymmetric, with the bottom side being the axis of revolution. At the top side, a slip boundary condition is imposed. The mean flow is quiescent with constant and unitary speed of sound and mean density. The frequency is chosen to be He = 2π, selected so that one complete period is contained within the computational domain.
The L1-norm of the error for the acoustic pressure is evaluated as where N dof is the number of degrees of freedom for each configuration. The exact solution to this problem isp i,exact = exp (−i He x) . Figure C.16 shows that for the less resolved configuration presented in this paper (10 elements per wavelength with second-order polynomials, p = 2) the error per wavelength is O(10 −3 ).