1. Introduction and objectives
Noise generated inside combustion chambers is receiving increasing attention from the scientific community as its relative contribution to global aircraft noise is growing. Waves generated inside the combustion chamber have to propagate through the turbine stages to contribute to the global noise at the outlet. This results in two mechanisms of noise generation (Marble & Candel Reference Marble and Candel1977): direct noise, where the acoustic fluctuations inside the combustion chamber propagate through turbine stages to the outlet, and indirect noise, which is generated by the acceleration of entropy waves (or hot spots) and vorticity waves through the turbine stages. Direct noise was initially addressed in the work of Strahle (Reference Strahle1972) and Hassan (Reference Hassan1974), and more recently by Rajaram & Lieuwen (Reference Rajaram and Lieuwen2009), Swaminathan et al. (Reference Swaminathan, Xu, Dowling and Balachandran2011), Talei, Brear & Hawkes (Reference Talei, Brear and Hawkes2011) and Ihme & Pitsch (Reference Ihme and Pitsch2012); indirect noise has been studied by Ffowcs Williams & Howe (Reference Ffowcs Williams and Howe1975), Howe (Reference Howe1975), Sinai (Reference Sinai1980) and Howe (Reference Howe2010). Since the noise generated has to propagate through the turbine stages, the simulation of the propagation of all these waves (acoustic, entropy and vorticity) through non-uniform flows is essential to correctly predict both mechanisms of noise generation.
The propagation of waves through the non-uniform flow downstream of the combustion chamber also generates reflected acoustic waves. These waves propagate upstream and interact with the flame, leading to an unsteady heat release which induces more acoustic and entropy waves. This can lead to a positive feedback, giving rise to thermoacoustic instabilities (Rayleigh Reference Rayleigh1894; Crocco Reference Crocco1969). If the acoustic energy generated by the flame–acoustics coupling is larger than the acoustic losses (through inflow and outflow boundaries, acoustic dissipation, etc.), the acoustic energy will increase, leading to an unstable thermoacoustic mode. The boundaries of the combustion chamber play a major role in the stability of the thermoacoustic mode: on the one hand, acoustic energy is lost through them, while on the other, acoustic energy can be generated at the boundaries when an entropy wave goes through them, contributing to the tendency for instability (Goh & Morgans Reference Goh and Morgans2013). For these reasons, the boundaries must be modelled accurately.
To compute the propagation of waves through a quasi-one-dimensional nozzle downstream of the combustion chamber, Marble & Candel (Reference Marble and Candel1977) assumed the nozzle to be short compared to the wavelength of the perturbations. This is known as the ‘compact assumption’ and provides both the reflection and the transmission coefficients of acoustic and entropy waves propagating through the nozzle in the limit of low frequencies. This limit has been found to be restrictive for non-zero frequencies (Leyko, Nicoud & Poinsot Reference Leyko, Nicoud and Poinsot2009), and different approaches have been developed to circumvent this limitation. One approach deals with a specified nozzle geometry; this was done by Marble & Candel (Reference Marble and Candel1977), who considered a linear mean velocity profile in a one-dimensional nozzle. This approach was then extended by Moase, Brear & Manzie (Reference Moase, Brear and Manzie2007) and Giauque, Huet & Clero (Reference Giauque, Huet and Clero2012), who considered a piecewise-linear mean velocity profile. A second approach was adopted by (Stow, Dowling & Hynes Reference Stow, Dowling and Hynes2002), who performed an asymptotic expansion in terms of frequency of the linearised Euler equations (LEE) of a choked nozzle. This yields an equivalent nozzle length, allowing a first-order correction to the phase for the reflection coefficient of choked nozzles. The method was implemented in a 2-D configuration, allowing for the propagation of 2-D circumferential modes. This ‘equivalent nozzle length’ method was further developed by Goh & Morgans (Reference Goh and Morgans2011), who obtained the phase correction of the transmission coefficient of a choked nozzle (with and without a shock wave) for both acoustic and entropy waves in the case of the one-dimensional plane mode.
The ‘equivalent nozzle length’ method considers only choked nozzles, and corrects only the phase of the reflection and transmission coefficients up to first order in frequency, since the modulus is computed using the ‘compact’ low-frequency assumption. A more complete approach was taken by Duran & Moreau (Reference Duran and Moreau2013), in which the LEE were solved in the one-dimensional case using the Magnus expansion (Magnus Reference Magnus1954; Blanes et al. Reference Blanes, Casas, Otero and Ros2009). This method provides the transfer functions of any quasi-one-dimensional geometry and the flow condition up to any frequency, correcting both the phase and modulus, for one-dimensional plane waves. The objective of this paper is to extend the Magnus expansion method to a two-dimensional annular configuration including circumferential waves. This is critical as aero-engine combustion chambers are typically annular. At the same time, this extension will allow the inclusion of vorticity waves in the model. These waves have been shown to generate indirect noise (Kings & Bake Reference Kings and Bake2010) and to have an influence on combustion instabilities (Li & Sun Reference Li and Sun2014).
Sketch of the annular configuration studied.
The rest of this article is organised as follows: § 2 presents the mathematical model to solve for the propagation of waves through the annular nozzle; § 3 briefly describes two numerical methods used to make comparisons with the mathematical model, and § 4 presents some of the results obtained. Conclusions are drawn in § 5.
2. Mathematical model
The aim of this section is to extend the one-dimensional solution of the LEE obtained by Duran & Moreau (Reference Duran and Moreau2013) using the Magnus expansion (Magnus Reference Magnus1954; Blanes et al.
Reference Blanes, Casas, Otero and Ros2009) to the more general case of an annular thin nozzle in a 2-D configuration, in which circumferential modes are considered. The annular duct is shown in figure 1, where cylindrical coordinates (
$x$
,
$r$
and
${\it\theta}$
) are used. For modelling purposes, the duct is considered to be thin, and therefore the mean radius will be used as the radius of the thin gap (
$R=(R_{out}(x)+R_{in}(x))/2$
) and the gap height,
$h(x)=R_{out}(x)-R_{in}(x)$
, will be considered small (
$h/R\ll 1$
). This will be called the ‘narrow-gap assumption’, and it has a double effect. On the one hand, the radial modes of the annular duct will be neglected, limiting the solution to frequencies lower than the first radial cut-off frequency (which is generally large as the ratio
$h/R$
is small for typical aero-engines). This limit is considered valid in the context of combustion noise and combustion instabilities, since frequencies are rarely high enough to contain radial modes in typical annular combustion chambers. On the other hand, the narrow-gap assumption will also have an effect on the mean flow, where radial effects will also be neglected. These effects can be particularly important in the case of convective waves. The validity of this assumption on the mean flow will be partially assessed using a numerical simulation (full-annular code; see § 3). It will also be assumed that the mean radius
$R$
is constant. This assumption is made to simplify the mathematical formulation, in line with the work of Stow et al. (Reference Stow, Dowling and Hynes2002), even though the mean radius of typical turbine stages decreases with
$x$
in general. The Euler equations integrated over
$r$
read
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left({\it\rho}h\right)+\frac{\partial }{\partial x}\left({\it\rho}uh\right)+\frac{\partial }{R\partial {\it\theta}}\left({\it\rho}vh\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left(u\right)+u\frac{\partial }{\partial x}\left(u\right)+v\frac{\partial }{R\partial {\it\theta}}\left(u\right)+\frac{1}{{\it\rho}}\frac{\partial }{\partial x}\left(p\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left(v\right)+u\frac{\partial }{\partial x}\left(v\right)+v\frac{\partial }{R\partial {\it\theta}}\left(v\right)+\frac{1}{{\it\rho}}\frac{\partial }{R\partial {\it\theta}}\left(p\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left(s\right)+u\frac{\partial }{\partial x}\left(s\right)+v\frac{\partial }{R\partial {\it\theta}}\left(s\right)=0, & \displaystyle\end{eqnarray}$$
${\it\rho}$
,
$s$
and
$p$
are the thermodynamic variables, i.e. density, entropy and pressure, and
$\boldsymbol{u}=[u,v]$
is the velocity vector with axial and circumferential components (the radial component is neglected).
Neglecting the temporal derivatives gives the laminar steady 2-D flow through the annular duct,
$$\begin{eqnarray}\frac{\partial }{\partial x}\left(\overline{{\it\rho}}\,\overline{u}h\right)=0,\quad \overline{u}\frac{\partial }{\partial x}\left(\overline{u}\right)+\frac{1}{\overline{{\it\rho}}}\frac{\partial }{\partial x}\left(\overline{p}\right)=0,\quad \overline{u}\frac{\partial }{\partial x}\left(\overline{v}\right)=0,\quad \overline{u}\frac{\partial }{\partial x}\left(\overline{s}\right)=0,\end{eqnarray}$$
where
$\bar{\left(\phantom{0}\right)}$
represents the steady mean flow given by the Euler equations (and does not include effects such as turbulence). In this work, as opposed to Stow et al. (Reference Stow, Dowling and Hynes2002), a constant circumferential component of the mean velocity will be considered.
To solve for the propagation of acoustic, entropy and vorticity waves through the annular nozzle, the Euler equations (2.1) are linearised. To do so, the primitive variables are written as
$$\begin{eqnarray}{\it\rho}=\overline{{\it\rho}}\left(1+\hat{{\it\rho}}\right),\quad p=\overline{p}\left(1+{\it\gamma}\hat{p}\right),\quad u=\overline{u}\left(1+\hat{u} \right),\quad v=\overline{v}+c_{t}\hat{v},\quad s=\overline{s}+c_{p}{\hat{s}},\end{eqnarray}$$
where
$\hat{\left(\phantom{0}\right)}$
represents small dimensionless perturbations and
${\it\gamma}$
is the ratio of specific heats. The perturbation in
$v$
is reduced using
$c_{t}$
, the stagnation sound speed (
$c_{t}=\sqrt{{\it\gamma}rT_{t}}$
, where
$r$
is the gas constant), and not the mean velocity
$\overline{v}$
, to avoid singularities when
$\overline{v}=0$
. Substituting (2.3) into (2.1) and neglecting second-order terms gives
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left(\hat{p}\right)+u\frac{\partial }{\partial x}\left(\hat{p}+\hat{u} \right)+v\frac{\partial }{R\partial {\it\theta}}\left(\hat{p}\right)+c_{t}\frac{\partial }{R\partial {\it\theta}}\left(\hat{v}\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left(\hat{u} \right)+u\frac{\partial }{\partial x}\left(\hat{u} \right)+v\frac{\partial }{R\partial {\it\theta}}\left(\hat{u} \right)+\frac{c^{2}}{u}\frac{\partial }{\partial x}\left(\hat{p}\right)+\left(2\hat{u} -({\it\gamma}-1)\hat{p}-{\hat{s}}\right)\frac{\text{d}u}{\text{d}x}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left(\hat{v}\right)+u\frac{\partial }{\partial x}\left(\hat{v}\right)+v\frac{\partial }{R\partial {\it\theta}}\left(\hat{v}\right)+\frac{c^{2}}{c_{t}}\frac{\partial }{R\partial {\it\theta}}\left(\hat{p}\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left({\hat{s}}\right)+u\frac{\partial }{\partial x}\left({\hat{s}}\right)+v\frac{\partial }{R\partial {\it\theta}}\left({\hat{s}}\right)=0, & \displaystyle\end{eqnarray}$$
${\hat{s}}=\hat{p}-\hat{{\it\rho}}$
has been used to eliminate
$\hat{{\it\rho}}$
from the equations. In the LEE above, the
$\bar{\left(\phantom{0}\right)}$
symbol for the mean flow has been dropped for simplicity.
To use the Magnus expansion, the equations are recast using the fluctuating mass, total temperature, entropy and azimuthal velocity,
$I_{A}=\hat{{\dot{m}}}$
,
$I_{B}=\hat{T_{t}}$
,
$I_{C}={\hat{s}}$
and
$I_{D}=\hat{v}$
, which, when written as functions of the primitive variables, read
$$\begin{eqnarray}\displaystyle & I_{A}=\hat{p}+\hat{u} -{\hat{s}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle I_{B}=({\it\gamma}-1)\frac{M_{x}^{2}\hat{u} +\hat{p}+\displaystyle \frac{{\hat{s}}}{{\it\gamma}-1}}{\displaystyle 1+\frac{{\it\gamma}-1}{2}M^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & I_{C}={\hat{s}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & I_{D}=\hat{v}. & \displaystyle\end{eqnarray}$$
$M_{x}=u/c$
and
$M_{{\it\theta}}=v/c$
are the axial and azimuthal Mach numbers (with
$c=\sqrt{{\it\gamma}rT}$
the speed of sound) and
$M=\sqrt{M_{x}^{2}+M_{{\it\theta}}^{2}}$
is the total Mach number. It is useful to first recast the mean axial-velocity gradient in terms of the Mach number gradient. This is done using the gradient of the square of the mean velocity modulus (
$w$
), 
$$\begin{eqnarray}\frac{\text{d}w^{2}}{\text{d}x}=\frac{\text{d}u^{2}}{\text{d}x}+\frac{\text{d}v^{2}}{\text{d}x},\end{eqnarray}$$
and, knowing that
$v$
is constant,
$$\begin{eqnarray}u\frac{\text{d}u}{\text{d}x}=Mc\frac{\text{d}(Mc)}{\text{d}x}.\end{eqnarray}$$
Using
$$\begin{eqnarray}\frac{\text{d}c}{\text{d}x}=-\frac{1}{2}\frac{({\it\gamma}-1)Mc}{\displaystyle 1+\frac{{\it\gamma}-1}{2}M^{2}}\frac{\text{d}M}{\text{d}x},\end{eqnarray}$$
the gradient of
$u$
reads
$$\begin{eqnarray}\frac{\text{d}u}{\text{d}x}=\frac{M}{M_{x}}\frac{c}{\displaystyle 1+\frac{{\it\gamma}-1}{2}M^{2}}\frac{\text{d}M}{\text{d}x}.\end{eqnarray}$$
The gradients of the axial and azimuthal Mach numbers (
$M_{x}$
and
$M_{{\it\theta}}$
) are also obtained as functions of the absolute Mach number gradient, as they will be used to simplify the equations. They read
$$\begin{eqnarray}\frac{\text{d}M_{x}}{\text{d}x}=\frac{M}{M_{x}}\frac{\displaystyle 1+\frac{{\it\gamma}-1}{2}M_{x}^{2}}{\displaystyle 1+\frac{{\it\gamma}-1}{2}M^{2}}\frac{\text{d}M}{\text{d}x},\quad \frac{\text{d}M_{{\it\theta}}}{\text{d}x}=\frac{\displaystyle \frac{{\it\gamma}-1}{2}MM_{{\it\theta}}}{\displaystyle 1+\frac{{\it\gamma}-1}{2}M^{2}}\frac{\text{d}M}{\text{d}x}.\end{eqnarray}$$
Equations (2.4) are finally written in terms of (2.5) using the relations in (2.9) and (2.10), giving
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial I_{A}}{\partial t}+u\frac{\partial I_{A}}{\partial x}+v\frac{\partial I_{A}}{R\partial {\it\theta}}+c_{t}\frac{\partial I_{D}}{R\partial {\it\theta}}+\frac{u}{({\it\gamma}-1)M_{x}^{2}}\left[{\it\zeta}\frac{\partial I_{B}}{\partial x}-\frac{\partial I_{C}}{\partial x}\right]=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial I_{B}}{\partial t}+u\frac{\partial I_{B}}{\partial x}+v\frac{\partial I_{B}}{R\partial {\it\theta}}+c_{t}\frac{{\it\gamma}-1}{{\it\zeta}}\frac{\partial I_{D}}{R\partial {\it\theta}}+u\frac{{\it\gamma}-1}{{\it\zeta}}\left[\frac{\partial I_{A}}{\partial x}+\frac{\partial I_{C}}{\partial x}\right]=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial I_{C}}{\partial t}+u\frac{\partial I_{C}}{\partial x}+v\frac{\partial I_{C}}{R\partial {\it\theta}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial I_{D}}{\partial t}+u\frac{\partial I_{D}}{\partial x}+v\frac{\partial I_{D}}{R\partial {\it\theta}}+c_{t}\frac{M_{x}^{2}}{\left(M_{x}^{2}-1\right){\it\zeta}}\frac{\partial I_{A}}{R\partial {\it\theta}}-\frac{c_{t}}{({\it\gamma}-1)\left(M_{x}^{2}-1\right)}\frac{\partial I_{B}}{R\partial {\it\theta}} & \displaystyle \nonumber\\ \displaystyle & \displaystyle \quad +c_{t}\frac{M_{x}^{2}+1/({\it\gamma}-1)}{\left(M_{x}^{2}-1\right){\it\zeta}}\frac{\partial I_{C}}{R\partial {\it\theta}}=0, & \displaystyle\end{eqnarray}$$
${\it\zeta}=1+({\it\gamma}-1)M^{2}/2$
for simplicity. In matrix form, (2.11) reads 
$$\begin{eqnarray}\frac{\partial }{\partial t}\boldsymbol{I}+u\unicode[STIX]{x1D640}_{\unicode[STIX]{x1D639}}\frac{\partial }{\partial x}\boldsymbol{I}+[v\unicode[STIX]{x1D644}+c_{t}\unicode[STIX]{x1D640}_{{\it\theta}}]\frac{\partial }{R\partial {\it\theta}}\boldsymbol{I}=0\end{eqnarray}$$
where
$\boldsymbol{I}$
is the vector of unknowns,
$\boldsymbol{I}=[I_{A}~I_{B}~I_{C}~I_{D}]^{\text{T}}$
,
$\unicode[STIX]{x1D644}$
is the identity matrix and
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D640}_{\unicode[STIX]{x1D639}}(x)=\left[\begin{array}{@{}cccc@{}}1 & \displaystyle \frac{{\it\zeta}}{({\it\gamma}-1)M_{x}^{2}} & -\displaystyle \frac{1}{({\it\gamma}-1)M_{x}^{2}} & 0\\ \displaystyle \frac{{\it\gamma}-1}{{\it\zeta}} & 1 & \displaystyle \frac{{\it\gamma}-1}{{\it\zeta}} & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\end{array}\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D640}_{{\it\theta}}(x)=\left[\begin{array}{@{}cccc@{}}0 & 0 & 0 & 1\\ 0 & 0 & 0 & \displaystyle \frac{{\it\gamma}-1}{{\it\zeta}}\\ 0 & 0 & 0 & 0\\ \displaystyle \frac{M_{x}^{2}}{{\it\zeta}\left(M_{x}^{2}-1\right)} & \displaystyle \frac{-1}{({\it\gamma}-1)\left(M_{x}^{2}-1\right)} & \displaystyle \frac{M_{x}^{2}+1/({\it\gamma}-1)}{{\it\zeta}\left(M_{x}^{2}-1\right)} & 0\end{array}\right]. & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D640}_{\unicode[STIX]{x1D639}}$
and
$\unicode[STIX]{x1D640}_{{\it\theta}}$
are functions of the
$x$
coordinate due to the axial dependence of the Mach number. Let us now consider harmonic waves of the form
$\exp \left(\text{i}wt+\text{i}n{\it\theta}\right)$
. Using the dimensionless numbers
${\it\Omega}=wL_{n}/c_{t}$
and
${\it\Omega}_{c}=nL_{n}/R$
, the dimensionless space coordinate
${\it\xi}=x/L_{n}$
(where
$L_{n}$
is the axial length of the nozzle) and the dimensionless velocity components
$U=u/c_{t}$
and
$V=v/c_{t}$
, the equation can be expressed as 
$$\begin{eqnarray}\frac{\partial }{\partial {\it\xi}}\boldsymbol{I}=\unicode[STIX]{x1D63C}\boldsymbol{I},\quad \text{with}~\unicode[STIX]{x1D63C}=\frac{-\text{i}}{U}[\unicode[STIX]{x1D640}_{\unicode[STIX]{x1D639}}]^{-1}[{\it\Omega}\unicode[STIX]{x1D644}+{\it\Omega}_{c}\left(\unicode[STIX]{x1D640}_{{\it\theta}}+V\unicode[STIX]{x1D644}\right)],\end{eqnarray}$$
which is equivalent to the equation obtained by Duran & Moreau (Reference Duran and Moreau2013) but with a fourth equation in
$I_{D}\equiv \hat{v}$
and an extra term in the matrix
$\unicode[STIX]{x1D63C}$
to account for circumferential waves (proportional to
${\it\Omega}_{c}$
).
Equation (2.15) can be solved using the Magnus expansion (Magnus Reference Magnus1954; Blanes et al. Reference Blanes, Casas, Otero and Ros2009), as shown for the one-dimensional case. The method considers a solution of the form
$$\begin{eqnarray}\boldsymbol{I}({\it\xi})=[\exp [\unicode[STIX]{x1D63D}({\it\xi})]]\boldsymbol{I}_{0}\quad \text{with}~\unicode[STIX]{x1D63D}({\it\xi})=\mathop{\sum }_{k=1}^{\infty }\unicode[STIX]{x1D63D}^{(k)}({\it\xi}),\end{eqnarray}$$
where
$\unicode[STIX]{x1D63D}^{(k)}({\it\xi})$
are the terms of the Magnus expansion (each of order
$O({\it\Omega},{\it\Omega}_{c})^{k}$
) and are obtained through a recursive procedure as explained by Blanes et al. (Reference Blanes, Casas, Otero and Ros2009). For convergence purposes, the series may have to be computed separately for different sections of the nozzle and then all the results multiplied together as explained in Blanes et al. (Reference Blanes, Casas, Otero and Ros2009) and Duran & Moreau (Reference Duran and Moreau2013).
Once the solution for
$\boldsymbol{I}({\it\xi})$
is found, its value at the inlet (
$\boldsymbol{I}_{0}$
) is prescribed using the boundary conditions to impose the incoming forcing waves. The method followed is similar to that used by Duran & Moreau (Reference Duran and Moreau2013), but some differences arise due to the presence of circumferential waves. The first step is to rearrange
$\boldsymbol{I}$
into incoming (known) and outgoing (unknown) waves. The relation can be written as
$\boldsymbol{I}=\unicode[STIX]{x1D63F}_{W}^{I}\boldsymbol{W}$
, where
$\boldsymbol{W}=[w^{+}~w^{-}~w^{v}~w^{s}]^{\text{T}}$
is a vector containing the two acoustic waves, an entropy wave and a vorticity wave. These waves are eigenvectors of the LEE (2.4), and their eigenvalues are related to their propagation speeds. In this way, the eigenvalue will determine whether a given wave (eigenvector) is incoming, and therefore has to be imposed at the boundary, or outgoing. In typical constant-Mach-number flows, the four different waves propagate independently and are said to be uncoupled. In a flow with non-homogeneous Mach number, such as the one considered here, the above is not true, and the propagation of entropy and vorticity waves is said to be coupled to the acoustic waves, since they generate indirect noise. To obtain matrix
$\unicode[STIX]{x1D63F}_{W}^{I}$
(and hence the relation between
$\boldsymbol{I}$
and
$\boldsymbol{W}$
), the primitive variables
$\boldsymbol{P}=[\hat{p}~\hat{u} ~\hat{v}~{\hat{s}}]^{\text{T}}$
are written as
$\boldsymbol{P}=\tilde{\boldsymbol{P}}\exp [\text{i}wt+\text{i}k_{x}x+\text{i}n{\it\theta}]$
and substituted into (2.4) to give an eigenvalue problem. The eigenvalues give the axial components of the dimensionless wavevector
$\tilde{k_{x}}=k_{x}c/w$
, which are obtained in terms of the mean flow and
${\it\mu}=(nc)/(Rw)$
, the ratio of the axial to azimuthal wavelengths (
${\it\mu}={\it\Omega}_{c}/({\it\Omega}\sqrt{{\it\zeta}})$
). For hydrodynamic waves (entropy and vorticity), the wavevector reads
$$\begin{eqnarray}\tilde{k_{x}}^{s,v}=-\frac{1+M_{{\it\theta}}{\it\mu}}{M_{x}}.\end{eqnarray}$$
For the acoustic waves, the wavevector is obtained through
$$\begin{eqnarray}\tilde{k_{x}}^{\pm }=\frac{M_{x}\left(1+M_{{\it\theta}}{\it\mu}\right)}{1-M_{x}^{2}}\mp \frac{{\it\Delta}^{1/2}}{1-M_{x}^{2}},\quad \text{where}~{\it\Delta}=[1+M_{{\it\theta}}{\it\mu}]^{2}-[1-M_{x}^{2}]{\it\mu}^{2}.\end{eqnarray}$$
When
${\it\Delta}<0$
, the axial component of the wavevector is complex, representing decaying waves. In this case
$\pm$
is used instead (waves decay in the propagation direction). The decay rate of an evanescent wave (in the absence of any other effect) is given by
$$\begin{eqnarray}\exp [\text{Im}(k_{x}^{+})x]=\exp \left[-\frac{\left[\left(1-M_{x}^{2}\right){\it\Omega}_{c}^{2}-\left({\it\Omega}\sqrt{{\it\zeta}}+M_{{\it\theta}}{\it\Omega}_{c}\right)^{2}\right]^{1/2}}{1-M_{x}^{2}}{\it\xi}\right].\end{eqnarray}$$
The eigenvectors of the system give a relation between the primitive variables and the waves,
$\boldsymbol{P}=\unicode[STIX]{x1D63F}_{W}^{P}\boldsymbol{W}$
, where
$$\begin{eqnarray}\unicode[STIX]{x1D63F}_{W}^{P}=\left(\begin{array}{@{}cccc@{}}1 & 1 & 0 & 0\\ {\it\eta}^{+} & {\it\eta}^{-} & {\it\Omega}_{c} & 0\\ {\it\beta}^{+} & {\it\beta}^{-} & -M_{x}{\it\Omega}\tilde{k_{x}}^{v} & 0\\ 0 & 0 & 0 & 1\end{array}\right),\quad \text{with }\begin{array}{@{}c@{}}\displaystyle {\it\eta}^{\pm }=-\frac{\tilde{k_{x}}^{\pm }}{M_{x}\left(1+M_{x}\tilde{k_{x}}^{\pm }+M_{{\it\theta}}{\it\mu}\right)},\\ \displaystyle {\it\beta}^{\pm }=-\frac{{\it\mu}}{\sqrt{{\it\zeta}}\left(1+M_{x}\tilde{k_{x}}^{+}+M_{{\it\theta}}{\it\mu}\right)}.\end{array}\end{eqnarray}$$
Note that, since the eigenvectors of any matrix are determined only up to their direction, the relation between
$\boldsymbol{P}$
and
$\boldsymbol{W}$
is not unique, as it depends on the magnitude considered. In this work the relation given by matrix
$\unicode[STIX]{x1D63F}_{W}^{P}$
of (2.20) will be used to define
$\boldsymbol{W}$
. Finally, using (2.5), matrix
$\unicode[STIX]{x1D63F}_{W}^{I}$
(relating
$\boldsymbol{I}$
to
$\boldsymbol{W}$
) reads
$$\begin{eqnarray}\unicode[STIX]{x1D63F}_{W}^{I}=\left(\begin{array}{@{}cccc@{}}1+{\it\eta}^{+} & 1+{\it\eta}^{-} & {\it\Omega}_{c} & -1\\ \displaystyle \frac{{\it\gamma}-1}{{\it\zeta}}\left(1+M_{x}^{2}{\it\eta}^{+}\right) & \displaystyle \frac{{\it\gamma}-1}{{\it\zeta}}\left(1+M_{x}^{2}{\it\eta}^{-}\right) & \displaystyle \frac{{\it\gamma}-1}{{\it\zeta}}M_{x}^{2}{\it\Omega}_{c} & \displaystyle \frac{1}{{\it\zeta}}\\ 0 & 0 & 0 & 1\\ {\it\beta}^{+} & {\it\beta}^{-} & -M_{x}{\it\Omega}\tilde{k_{x}}^{v} & 0\end{array}\right).\end{eqnarray}$$
The matrix
$\unicode[STIX]{x1D63F}_{W}^{I}$
is now used to impose the forcing waves in order to obtain
$\boldsymbol{I}_{0}$
, as explained by Duran & Moreau (Reference Duran and Moreau2013). It should be noticed that for the case where
${\it\Omega}_{c}=0$
and
$M_{{\it\theta}}=0$
are considered, the solution is identical to that in the one-dimensional case. Once obtained, the waves at any point (particularly those at the inlet and outlet of the nozzle) can be calculated using the Magnus expansion (2.16) to obtain the vector
$\boldsymbol{I}$
at the desired location and the matrix
$\unicode[STIX]{x1D63F}_{W}^{I}$
to transform the solution into the vector
$\boldsymbol{W}$
.
3. Numerical methods
Two numerical methods will be used to perform comparisons with the analytical solution. The first numerical method, hereafter referred to as narrow-gap, solves (2.4) for a given mean flow in a one-dimensional grid and considers an angular dependence of the perturbations in the form
$\exp \left(\text{i}n{\it\theta}\right)$
. Since the equations (2.4) are being solved, the narrow-gap assumption is used and no radial modes are considered in this numerical method. The code uses an eighth-order centred numerical scheme with 4000 grid points and high-order explicit filtering to avoid numerical fluctuations. Time-harmonic disturbances are assumed, and the solution is advanced using pseudo-time-stepping and a fourth-order Runge–Kutta method. The solution has been shown to be grid-independent.
The second code, hereafter referred to as full-annular, was used by Stow et al. (Reference Stow, Dowling and Hynes2002), where a detailed explanation of the method can be found. The objective of this code is to evaluate the influence of the narrow-gap assumption made in both the narrow-gap code and the Magnus expansion method. The geometry considered is therefore the full annular nozzle of figure 1, solved while considering azimuthal symmetry. The grid is 2-D and the flow is solved using
$x$
and
$r$
coordinates with azimuthal symmetry. The 2-D Euler equations are first solved in the annular nozzle with no azimuthal velocity to obtain the mean flow as a function of the axial and radial coordinates, hence without the narrow-gap assumption. This is done using a second-order finite volume code as explained in Stow et al. (Reference Stow, Dowling and Hynes2002), with 160 grid points in the axial direction and 80 in the radial direction. The mean flow is then used by a linearised code (running on the same grid and using the same numerical methods) to solve for the propagation of waves, assuming an angular dependence of the perturbations in the form
$\exp \left(\text{i}n{\it\theta}\right)$
and using pseudo-time-stepping for each frequency. This method does not impose the narrow-gap assumption, as variations in the radial direction are considered. The three methods (the analytical Magnus expansion and the two numerical methods narrow-gap and full-annular) will be compared.
4. Transfer functions through annular nozzles
In this section the transfer functions of an annular nozzle will be presented. Both the reflection and the transmission coefficients will be used, which are defined, respectively, as the reflected and transmitted acoustic wave divided by the incoming perturbation. For the acoustic transfer function, for instance, the reflection coefficient is
$w_{1}^{-}/w_{1}^{+}$
while the transmission coefficient is
$w_{2}^{+}/w_{1}^{+}$
(where subscript 1 refers to the inlet and subscript 2 to the outlet); similarly for the entropy wave, with the incoming perturbation
$w_{1}^{+}$
replaced by
$w_{1}^{s}$
. The two acoustic waves (
$w^{+}$
and
$w^{-}$
) and the entropy and vorticity waves (
$w^{s}$
and
$w^{v}$
) are related to the primitive variables through matrix (2.20) and to vector
$\boldsymbol{I}$
through matrix (2.21). These transfer functions can be viewed as the response of the nozzle to a ‘unitary’ perturbation, since the problem is linear. The unitary acoustic wave at the inlet, for example, will imply
$w_{1}^{+}=1$
as an inlet forcing wave. However, it should be noted that the definition of the waves is linked to (2.20), which is the eigenvector matrix. Since eigenvectors are defined only up to their direction (as explained previously), the definition of the waves is not universal and may vary from one study to another; therefore the transfer functions discussed here may differ according to the definition used for the acoustic, entropy and vorticity waves. In this work (2.20) is used throughout.
The purpose of this section is to compare the analytical method with the numerical simulations and show that the Magnus expansion method is able to predict the propagation of waves through annular nozzles for non-compact frequencies. One representative geometry will be used for the study as an example, although the method can be applied to any annular nozzle geometry, for any flow configuration (subsonic, choked etc.), provided that the narrow-gap annular nozzle assumption is valid (
$h\ll R$
).
The present work will focus on results for the circumferential waves (
$n\neq 0$
), which is the novel contribution of this paper, as use of the Magnus expansion for plane modes has been discussed previously (Duran & Moreau Reference Duran and Moreau2013). The geometry considered is sketched in figure 1, while figure 2 shows the area of the nozzle as a function of the axial coordinate, normalised by the critical section
$A^{\ast }$
. The geometry is similar to that used by Stow et al. (Reference Stow, Dowling and Hynes2002), but with a divergent section that is long compared to the convergent section. This is typical of aero-engine designs, in order to avoid flow separation in the divergence. The total pressure in the nozzle is 2.16 bars and the total temperature is 986 K; the inlet and outlet axial Mach numbers are 0.29 and 2.2, respectively. Figure 3(a) shows the pressure field through a cross-section of the nozzle at constant
${\it\theta}$
(obtained using the Euler solver). The mean Mach number averaged over
$r$
(computed using the full-annular code) is plotted in figure 3(b). The mean Mach number is used for the Magnus expansion and the narrow-gap code, as these approaches assume a narrow-gap annular nozzle, whereas the full-annular code considers the radial dependence of the mean flow shown in figure 3(a). In §§ 4.1–4.3 the mean flow will be axial,
$M_{{\it\theta}}=0$
. The effect of non-zero azimuthal mean flow will be discussed in § 4.4.
Area of the nozzle, computed as
$A=h{\rm\pi}(R_{out}^{2}-R_{in}^{2})$
.
Mean flow through the annular nozzle, computed using the full-annular code: (a) cross-section of the annular nozzle showing the pressure field (in bars); (b) mean Mach number averaged over
$r$
.
4.1. Acoustic transfer functions
Figure 4 shows the acoustic transfer functions for
$n=1$
, equivalent to
${\it\Omega}_{c}=1$
. It can be seen that the analytical method produces results that are in perfect agreement with both numerical simulations: the analytical method using the Magnus expansion is able to compute both the reflection and the transmission coefficients of annular thin nozzles at a fraction of the cost of a numerical simulation. Results show a cut-off frequency near
${\it\Omega}<{\it\Omega}_{c}\left(1-M_{x}^{2}\right)^{1/2}\!/\sqrt{{\it\zeta}}$
(from
${\it\Delta}=0$
in (2.18)). For frequencies lower than the cut-off frequency, the acoustic waves are evanescent in the low-Mach-number inlet region. It is worth emphasising that acoustic waves are evanescent only in the inlet region of the nozzle, where the Mach number is low, and propagate normally for larger values of the Mach number, as can be seen from the discriminant of (2.18). This means that the waves decay only in a small inlet region of the nozzle, with a decay rate given by (2.19). This is shown in figure 5, where the absolute values of the waves propagating through the nozzle are represented for two different frequencies, one cut-off and one cut-on. The peak value seen close to
${\it\xi}=0.3$
for
${\it\Omega}=0.2$
corresponds to the value of the Mach number for which
${\it\Delta}=0$
, the point where the evanescent acoustic waves become cut-on. A singularity is found at this point, as the axial component of the wavevector is the same for both acoustic waves (as seen in (2.18)): one cannot distinguish between the two acoustic waves, as they both propagate in the same direction and at the same speed, and consequently
$\unicode[STIX]{x1D63F}_{W}^{P}$
is singular. Figure 5 also shows that the propagation of acoustic waves is uncoupled from the propagation of the entropy and vorticity waves.
Acoustic transfer functions of the nozzle for
$n=1$
: (a) reflection coefficient; (b) transmission coefficient.
Waves propagating through the nozzle with a unitary acoustic wave for
$n=1$
: (a)
${\it\Omega}=0.2$
; (b)
${\it\Omega}=1.5$
.
Figure 4(b) shows that the transmission coefficient of acoustic waves increases strongly at low frequencies, despite these frequencies being cut-off. It should be noted that, although the acoustic waves are evanescent, they do not decay by more than 10 %, as seen from (2.19) and shown in figure 5(a). The increase at low frequencies of the transmitted wave is due to the definition of the acoustic waves. These waves are defined to have
$\hat{p}\propto w^{\pm }$
and
$\hat{u} \propto {\it\eta}^{\pm }w^{\pm }$
, as seen from (2.20). For evanescent waves the absolute value of
${\it\eta}^{\pm }$
is large, meaning that a unitary acoustic wave (carrying a unitary pressure fluctuation) contains large velocity fluctuations.
The acoustic transfer functions for
$n=4$
(
${\it\Omega}_{c}=4$
) are shown in figure 6. Again, the analytical Magnus expansion agrees with both numerical simulations. Figure 7 shows the evolution through the nozzle of the absolute value of the waves for the same two frequencies as before, but with
$n=4$
. Both frequencies exhibit evanescent waves in the inlet region, and it can be seen that the decay is significantly stronger (by up to 45 % in the inlet region only) than in the previous example with
$n=1$
.
Acoustic transfer functions of the nozzle for
$n=4$
: (a) reflection coefficient; (b) transmission coefficient.
Waves propagating through the nozzle with a unitary acoustic wave for
$n=4$
: (a)
${\it\Omega}=0.2$
; (b)
${\it\Omega}=1.5$
.
4.2. Entropy transfer functions
Figure 8 shows the reflected and transmitted acoustic waves generated by a unitary entropy wave. The transmitted waves are responsible for the so-called indirect combustion noise (Marble & Candel Reference Marble and Candel1977; Bake et al. Reference Bake, Richter, Muhlbauer, Kings, Rohle, Thiele and Noll2009; Howe Reference Howe2010), while the reflected waves can contribute to thermoacoustic instabilities (Goh & Morgans Reference Goh and Morgans2013; Motheau, Nicoud & Poinsot Reference Motheau, Nicoud and Poinsot2014). Again, good agreement is found between the analytical Magnus expansion and the numerical methods. The cut-off frequency is also found to be a discontinuity of the solution, presenting strong reflected acoustic waves and low levels of entropy noise. It can be seen that the propagation of the entropy wave generates acoustic waves, demonstrating the one-way coupling between entropy waves and acoustic waves in non-homogeneous flows. This coupling is said to be one-way only, since acoustic waves do not generate entropy waves, as seen in § 4.1.
Entropy transfer functions of the nozzle for
$n=1$
: (a) reflection coefficient; (b) transmission coefficient.
Waves propagating through the nozzle with a unitary entropy wave for
$n=1$
: (a)
${\it\Omega}=0.2$
; (b)
${\it\Omega}=1.5$
.
The evolution of the absolute value of the waves through the nozzle is plotted in figure 9 for two different reduced frequencies. In contrast to the case of the acoustic wave, there is a generation of vorticity due to the presence of an entropy wave, suggesting that the propagation of the entropy wave is also coupled to that of the vorticity wave. The explanation of this phenomenon can be found when one writes an evolution equation for the vorticity wave. Writing
${\it\varsigma}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{v}$
for the vorticity, the
$x$
and
$y$
momentum equations of (2.1) combine to give
$$\begin{eqnarray}\frac{\text{D}{\it\varsigma}}{\text{D}t}+{\it\varsigma}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{R\partial {\it\theta}}\right)+\frac{1}{{\it\rho}^{2}}\left(\frac{\partial {\it\rho}}{\partial x}\frac{\partial p}{R\partial {\it\theta}}-\frac{\partial {\it\rho}}{R\partial {\it\theta}}\frac{\partial p}{\partial x}\right)=0.\end{eqnarray}$$
The equation for the fluctuating vorticity,
${\it\varsigma}^{\prime }$
, is obtained by substituting
${\it\varsigma}=\overline{{\it\varsigma}}+{\it\varsigma}^{\prime }$
into (4.1) to get
$$\begin{eqnarray}\frac{\text{D}{\it\varsigma}^{\prime }}{\text{D}t}+{\it\varsigma}^{\prime }\frac{\text{d}u}{\text{d}x}-\left[\frac{Mc^{2}}{{\it\zeta}}\frac{\text{d}M}{\text{d}x}\right]\frac{\partial {\hat{s}}}{R\partial {\it\theta}}=0,\end{eqnarray}$$
where again the
$\bar{\left(\phantom{0}\right)}$
symbol for the mean has been dropped for simplicity. Upon inverting (2.20), the dimensionless vorticity wave reads
$$\begin{eqnarray}w^{v}=\frac{\text{i}{\it\varsigma}^{\prime }L_{n}}{u\left({\it\mu}^{2}+(\tilde{k}_{x}^{v})^{2}\right)}.\end{eqnarray}$$
Equation (4.2) shows that vorticity is not simply convected, as two extra terms appear in the convection equation of
${\it\varsigma}^{\prime }$
. The first is due to the vortex stretching, and it is proportional to the mean flow gradient. The second comes from the baroclinic torque: the fluctuating density (due to the entropy wave) is misaligned with respect to the mean pressure gradient when
$n\neq 0$
; this causes a torque force, since the mean pressure gradient induces a stronger acceleration of the lower-density regions of the flow, resulting in a non-zero azimuthal gradient of
$\hat{u}$
, and consequently vorticity. This vorticity wave is relevant, as it can generate indirect vorticity noise when further propagated (see § 4.3). It is therefore interesting to plot the vorticity wave generated by an entropy wave at different frequencies. This is done in figure 10, showing that the vorticity wave generated is significant mostly at low frequencies.
Vorticity generated by an entropy wave: (a)
$n=1$
; (b)
$n=4$
.
Figure 11 plots the entropy transfer functions for
$n=4$
. The acoustic waves generated are, in this particular case, weaker than for the previous
$n=1$
case. A small discrepancy can be seen between the full-annular code and the other two methods at
${\it\Omega}/2{\rm\pi}\sim 0.6$
. This is probably due to the inclusion of the radial effects on the propagation of the entropy wave in the full-annular code, as they are not included in the Magnus expansion method or the narrow-gap code.
Entropy transfer functions of the nozzle for
$n=4$
: (a) reflection coefficient; (b) transmission coefficient.
4.3. Vorticity transfer functions
A different source of indirect noise is vorticity waves, as explained theoretically by Cumpsty & Marble (Reference Cumpsty and Marble1977) and shown experimentally by Kings & Bake (Reference Kings and Bake2010). These waves have also been found to affect combustion instabilities (Li & Sun Reference Li and Sun2014). Vorticity waves stretch due to the presence of the mean flow gradient, generating acoustic waves that travel upstream (reflected waves) or downstream (transmitted noise). It is therefore interesting to obtain the transfer functions of the nozzle when perturbed by a unitary vorticity wave. These transfer functions are plotted in figure 12, where the analytical method is shown to correctly predict the acoustic response of the nozzle computed numerically. The same transfer functions, but with
$n=4$
, are plotted in figure 13. In this case the acoustic waves generated are significantly larger, although this is only due to the normalisation criteria used for the vorticity wave (see (4.3)). Small levels of numerical dissipation can be seen for both convective waves (entropy and vorticity) at large values of
${\it\Omega}$
. The convective wavelength,
${\it\lambda}\sim cM_{x}/{\it\Omega}$
, is low due to the large values of
${\it\Omega}$
and the low Mach number at the inlet of the nozzle, and so small levels of numerical dissipation appear when a fixed number of grid points is used.
Vorticity transfer functions of the nozzle for
$n=1$
: (a) reflection coefficient; (b) transmission coefficient.
Vorticity transfer functions of the nozzle for
$n=4$
: (a) reflection coefficient; (b) transmission coefficient.
4.4. The effect of mean azimuthal velocity
Up to this point, the cases studied concerned a flow with zero mean azimuthal velocity; however, this is not always the case. The mean flow at the outlet of the combustion chamber is not perfectly axial, and may contain an azimuthal component. Turbine stages impose strong flow deviations inside the annular duct. In this section, the differences between the
$M_{{\it\theta}}=0$
and
$M_{{\it\theta}}\neq 0$
cases will be highlighted. To do so, we use the same geometry as above, but consider a non-zero inlet azimuthal velocity that is 30 % of the mean axial velocity (equivalent to a mean flow angle of
${\it\theta}_{0}=16.7^{\circ }$
at the inlet). The mean flow through the nozzle is then computed using (2.2) and the propagation of acoustic, entropy and vorticity waves simulated using both the Magnus expansion and the narrow-gap code.
Figure 14 shows the acoustic transfer functions of the nozzle for the
${\it\theta}_{0}=16.7^{\circ }$
case, compared with the previous
${\it\theta}_{0}=0^{\circ }$
case for
$n=4$
. Excellent agreement is found between the analytical method and the narrow-gap code, showing that the Magnus expansion is able to compute the propagation of waves through nozzles, even when a non-zero mean azimuthal velocity is present. Figure 15 shows the same comparison for the entropy transfer functions. It should be noticed that the main influence of non-zero mean azimuthal velocity is to change the cut-off frequency of the transfer function.
Acoustic transfer functions of the nozzle for
$n=4$
with
${\it\theta}_{0}=16.7^{\circ }$
: (a) reflection coefficient; (b) transmission coefficient.
Entropy transfer functions of the nozzle for
$n=4$
with
${\it\theta}_{0}=16.7^{\circ }$
: (a) reflection coefficient; (b) transmission coefficient.
5. Conclusions
An analytical method for solving for the one-dimensional propagation of waves, developed by Duran & Moreau (Reference Duran and Moreau2013) and based on the Magnus expansion, is extended in this article to a two-dimensional configuration with circumferential waves. Our approach includes four novel aspects not considered previously. First, it considers circumferential modes as well as plane modes; this is essential for typical annular combustors, in which such modes are also present. Second, it includes the effect of a non-zero mean flow angle, which is of significant relevance in turbine stages. Third, it includes the vorticity wave, which is responsible for part of the indirect combustion noise, as shown by Kings & Bake (Reference Kings and Bake2010). Finally, it shows that circumferential entropy waves generate vorticity when accelerated; this vorticity wave, shown to originate due to the baroclinic torque of the vorticity equation, is of first order and can contribute to indirect vorticity noise when further accelerated and decelerated, as seen in the literature and in § 4.3.
Our extension of the analytical method is able to accurately predict the transfer functions for acoustic, entropy and vorticity waves for any circumferential mode, as demonstrated by comparisons with numerical simulations. Although not shown in these results, the mathematical method presented here can also be applied to subsonic flows through nozzles and to choked flows with a shock wave (Stow et al. Reference Stow, Dowling and Hynes2002; Duran & Moreau Reference Duran and Moreau2013). Prediction of these transfer functions using simplified analytical methods is ideal for use in low-order modelling of both combustion instabilities (Stow & Dowling Reference Stow and Dowling2009; Poinsot & Veynante Reference Poinsot and Veynante2011; Li & Morgans Reference Li and Morgans2015) and combustion noise (Duran et al. Reference Duran, Moreau, Nicoud, Livebardon, Bouty and Poinsot2014).
Acknowledgement
This work was funded by the ERC starting grant ACOULOMODE (ref. 305410); this support is gratefully acknowledged.
Open access
