1. Introduction
Solutions of the unsteady response of ducts with area variations and sustaining a mean flow, i.e. nozzle flows, are of interest for a wide variety of industrial applications, including combustors, automotive exhausts, after-burners or supersonic air-intake diffusers. An early attempt to mathematically describe the transfer function of rocket nozzles was proposed by Tsien (Reference Tsien1952), motivated by emerging theoretical descriptions of combustion instabilities in rocket engines.
With the same motivation, Marble & Candel (Reference Marble and Candel1977) derived an analytical solution for the transfer function of both subcritical and supercritical nozzles excited by acoustic and/or entropy waves. This solution is valid in the compact (or zero-frequency) limit which assumes that the acoustic and entropy wavelengths significantly exceed the length scale over which the nozzle area change occurs. This compact solution was extended to non-zero frequencies by Stow, Dowling & Hynes (Reference Stow, Dowling and Hynes2002) and Goh & Morgans (Reference Goh and Morgans2011) using an asymptotic expansion of the linearised Euler equations in terms of frequency. Building on the same assumptions of Marble & Candel (Reference Marble and Candel1977), namely inviscid, isentropic and quasi-one-dimensional flow, Duran & Moreau (Reference Duran and Moreau2013) proposed a solution valid at any frequency based on the Magnus expansion method (Blanes et al. Reference Blanes, Casas, Oteo and Ros2009). This approach was later extended to account for incoming compositional inhomogeneities by Magri (Reference Magri2017), to annular nozzles by Duran & Morgans (Reference Duran and Morgans2015) and to multi-stream nozzles by Younes & Hickey (Reference Younes and Hickey2019). All the previous studies were developed in the context of indirect combustion noise (Ihme Reference Ihme2017), which is the noise generated when convective disturbances (entropy, vortical and compositional waves) are accelerated/decelerated. Indirect combustion noise is particularly relevant for high- Mach-number flows and low frequencies.
In parallel to the theoretical efforts of the combustion noise community, solutions for the acoustic field in nozzle flows were presented for low subsonic Mach numbers. Dokumaci (Reference Dokumaci1998) proposed a semi-analytical solution based on the WKB (Wentzel–Kramers–Brillouin) method which requires sufficiently high frequencies. Additionally, several analytical solutions of the acoustic field in nozzles with specific area profiles were obtained by Eisenberg & Kao (Reference Eisenberg and Kao1971), Easwaran & Munjal (Reference Easwaran and Munjal1992) and Subrahmanyam, Sujith & Lieuwen (Reference Subrahmanyam, Sujith and Lieuwen2001) using transformation of variables.
An assumption of all the aforementioned studies is the isentropicity of the flow. However, many practical applications require solutions for non-isentropic mean flows. For instance, De Domenico et al. (Reference De Domenico, Rolland, Rodrigues, Magri and Hochgreb2021) and Yang, Guzmán-Iñigo & Morgans (Reference Yang, Guzmán-Iñigo and Morgans2020) recently showed that, in the diverging portion of realistic nozzles, the flow can separate creating recirculation regions and, hence, turbulence both leading to mean flow non-isentropicity and strong deviations from isentropic acoustic theories. Another source of non-isentropicity can be heat transfer, as happens in combustors and heat exchangers. Yeddula, Gaudron & Morgans (Reference Yeddula, Gaudron and Morgans2021) showed that heat transfer can significantly affect the acoustics of nozzles. However, there is a lack of solutions to this problem, even in the compact limit. The only solution to date was proposed by Yeddula & Morgans (Reference Yeddula and Morgans2021) based on the WKB method and assumes high frequencies and low subsonic flow Mach numbers.
In this paper, we propose a general solution for the acoustic field produced by acoustic and/or entropy waves in a nozzle whose mean flow is non-isentropic due to steady heat transfer. The solution is based on the Magnus expansion method of Duran & Moreau (Reference Duran and Moreau2013) and is valid for both subcritical and supercritical flow conditions, and for any frequency within the limit of the quasi-one-dimensional assumption. The Magnus expansion when applied to nozzle flows has always so far had frequency as the expansion parameter. In this work, we include a measure of the flow non-isentropicity as an expansion parameter for the first time, allowing us to develop models for non-isentropic nozzle flows which are valid at all frequencies and Mach numbers. The non-isentropicity considered is due to axially varying steady heat transfer; the method could be extended to include other forms of non-isentropicity in future work.
This article is organised as follows. In § 2 the mathematical model is described and a solution is proposed based on the Magnus expansion method. The solution is applied to a converging–diverging nozzle in § 3 and conclusions are drawn in § 4.
2. Analysis
We consider a calorically perfect gas flowing through a nozzle as depicted in figure 1 with steady heat transfer which varies axially as 
$\bar {\dot {Q}}(x)$. The flow is taken to be inviscid and compressible, and a quasi-one-dimensional framework is adopted. Neglecting volumetric forces as well as thermal and mass diffusion, the conservation of mass, momentum, energy and equation of state can be written, respectively, as,
\begin{align} A \frac{ \partial \rho }{ \partial t } +\frac{\partial(\rho Au)}{\partial x} = 0, \quad \frac{ \partial u }{ \partial t } +u\frac{\partial{u}}{\partial x} + \frac{1}{\rho}\frac{\partial p}{\partial x} = 0, \quad \frac{ \partial {s} }{ \partial t } +{u}\frac{\partial{s}}{\partial x} = \frac{{R}_{g}}{p}{\bar{\dot{Q}}} \quad \text{and} \quad {p} = {\rho} {R}_{g} {T}, \end{align}
where 
${\rho }$ is density, 
${p}$ is pressure, 
${T}$ is temperature, 
${u}$ is axial velocity, 
${s}$ is entropy, 
${R}_{g}$ is the gas constant, 
$A$ is the cross-sectional area of the nozzle and 
${\bar {\dot {Q}}}$ denotes the steady volumetric heat source term.
Sketch of a nozzle of length 
$L$ exchanging heat with the surroundings. Here 
$A(x)$, 
$M(x)$, and 
$\bar {T}(x)$ represent the nozzle cross-sectional area, flow Mach number and mean temperature at any 
$x$ varying from 
${x}_{0}$ to 
${x}_{1}$.
We seek to retrieve the dynamics of unsteady, small-amplitude perturbations superimposed on a steady background mean flow. The linearisation principle allows the thermodynamic and flow variables to be decomposed into the sum of mean time-averaged and fluctuating time-dependant components, denoted by 
$\overline {(\ )}$ and 
${(\;)'}$, respectively, e.g 
$\rho = \bar {\rho } + \rho '$. These time-dependant quantities are further normalised as follows:
\begin{equation} \hat{\rho}=\dfrac{\rho'}{\bar{\rho}}, \quad \hat{p}=\dfrac{p'}{\gamma\bar{p}}, \quad \hat{u}=\dfrac{u'}{\bar{u}}, \quad \hat{s}=\left(\gamma - 1\right)\dfrac{s'}{\gamma R_g}, \end{equation}
with 
$\gamma$ denoting the adiabatic index. It is assumed that there are no heat fluctuations inside the nozzle (
${\dot {Q}}^{'} = 0$). The analysis could be extended to account for heat fluctuations, but would require an additional closure model linking these to flow fluctuations. For simplicity, this paper focuses on the primary effect of mean heat transfer.
By linearising the thermodynamic and flow variables in the conservation equations (2.1) and equating the mean quantities, we obtain
\begin{align}
\dfrac{1}{\bar{\rho}}\dfrac{\text{d}\bar{\rho}}{{\text{d}\kern0.06em x}}+ \dfrac{1}{\bar{u}}\dfrac{\text{d}\bar{u}}{\text{d}\kern0.06em x}
+ \dfrac{1}{A}\dfrac{\text{d}A}{\text{d}\kern0.06em x} = 0, \enspace
\bar{u}\dfrac{\text{d}\bar{u}}{{\text{d}\kern0.06em x}}
={-}\dfrac{1}{\bar{\rho}}\dfrac{\text{d}\bar{p}}{\text{d}\kern0.06em x},\enspace \dfrac{\gamma R_g}{\gamma -
1}\dfrac{\text{d}\bar{T}}{\text{d}\kern0.06em x} +
\bar{u}\dfrac{\text{d}\bar{u}}{{\text{d}\kern0.06em x}} =
\dfrac{\bar{\dot{Q}}}{\bar{\rho} \:\bar{u}} \enspace
\text{and} \enspace \bar{p} = \bar{\rho}R_g \bar{T}.
\end{align}
Similarly, the linearised forms of the Euler equations are obtained in the time domain by equating the first-order fluctuating quantities of the conservation equations (2.1) in their normalised form given by (2.2). They read
\begin{gather} \dfrac{ \partial \hat{p} }{ \partial t } + {\bar{u}} \dfrac{\partial}{\partial x}\left(\hat{p} + \hat{u} \right) ={-} \bar{\dot{Q}} \dfrac{\gamma - 1}{\gamma \bar{p}}\left(\hat{u} + \gamma\hat{p}\right), \end{gather}
\begin{gather} \dfrac{ \partial \hat{u}}{ \partial t } + {\bar{u}} \dfrac{\partial \hat{u}}{\partial x} + \dfrac{\bar{c}^2}{\bar{u}} \dfrac{\partial \hat{p}}{\partial x} + \dfrac{\text{d}{\bar{u}}}{\text{d}\kern0.06em x}\left[2\hat{u} - \left(\gamma-1\right)\hat{p} - \hat{s} \right] = 0, \end{gather}
\begin{gather} \dfrac{ \partial \hat{s}}{ \partial t } + {\bar{u}} \dfrac{\partial \hat{s}}{\partial x} ={-} \bar{\dot{Q}} \dfrac{\gamma - 1}{\gamma \bar{p}}\left(\hat{u} + \gamma\hat{p}\right), \end{gather}
with 
$\bar {c}$ denoting the adiabatic speed of sound. The above linearised Euler equations (LEE) (2.4)–(2.6) govern the perturbed flow in a nozzle sustaining a mean flow with steady heat transfer.
To employ the Magnus expansion, the linearised Euler equations need to be recast in terms of three variables that are invariants of the flow at zero frequency and with zero heat transfer. The normalised fluctuating mass flow rate (
$I_A = \hat {m}$), stagnation temperature (
$I_B = \hat {T}_t$) and entropy (
$I_C = \hat {s}$) are chosen as the three invariants, where,
\begin{equation} \hat{m} \equiv
\dfrac{m'}{\bar{m}}, \quad \hat{T}_t \equiv
\dfrac{T'_t}{\bar{T}_t},
\end{equation}
with 
$m = \rho A u$ and 
$T_t = T(1 + (({\gamma -1})/{2})({u}/{c})^2).$ In terms of the primitive variables 
$\hat {p}, \hat {u}$ and 
$\hat {s}$, these invariants read
\begin{equation} I_{A}=\hat{p}+\hat{u}-\hat{s}, \quad I_{B}=(\gamma-1) \dfrac{M^{2} \hat{u}+\hat{p}+\dfrac{\hat{s}}{\gamma-1}} {1+\dfrac{\gamma-1}{2} M^{2}}, \quad I_{C}=\hat{s}, \end{equation}
where 
$M=\bar {u}/\bar {c}$ is the local Mach number of the flow. By defining a vector of invariants such that 
$\boldsymbol{\mathsf{I}} = {[I_A \;I_B\; I_C]}^{\text {T}}$ and using 
$\zeta = {1+(({\gamma -1})/{2}) M^{2}},$ the linearised Euler equations can be recast in matrix form. The full derivation is outlined in Appendix A and the final equation reads
\begin{equation} \dfrac{\partial}{\partial t} \boldsymbol{\mathsf{I}} +\bar{u} \boldsymbol{\mathsf{E}}_{x} \frac{\partial}{\partial x} \boldsymbol{\mathsf{I}}+ \dfrac{\left(\gamma -1\right) \bar{\dot{Q}}}{2\zeta \gamma \bar{p}}\boldsymbol{\mathsf{E}}_{s} \boldsymbol{\mathsf{I}} = 0, \quad \text{where} \quad \boldsymbol{\mathsf{E}}_x=\left[\begin{array}{ccc} 1 & \dfrac{\zeta}{(\gamma-1) M^{2}} & -\dfrac{1}{(\gamma-1) M^{2}} \\ \dfrac{\gamma-1}{\zeta} & 1 & \dfrac{\gamma-1}{\zeta} \\ 0 & 0 & 1 \end{array}\right], \end{equation}
\begin{equation} \boldsymbol{\mathsf{E}}_s=\dfrac{1}{M^2 - 1}\left[\begin{array}{ccc} {M^2\left(1-\gamma\right)-2} & \zeta \dfrac{\gamma+1}{\gamma-1} & -\gamma M^2 - \dfrac{2\gamma}{\gamma-1} \\ 2\gamma\left(\gamma-1\right) + \dfrac{M^2\left(1-\gamma^2\right)}{\zeta} & -2\zeta\gamma + \dfrac{\gamma^2 M^2}{\gamma-1} & 2\gamma\left(\gamma-\dfrac{ M^2}{\zeta}\right) -\dfrac{\gamma\left(\gamma-1\right)}{\zeta}M^4\\ 2\zeta\left(\gamma M^2-1\right) & -2{\zeta}^2 & 2\zeta\gamma M^2 \end{array}\right]. \end{equation}
As observed in (2.9) and (2.10), the coefficients of the matrices are only functions of the axial mean flow and the mean rate of heat transfer per unit volume, 
$\bar {\dot {Q}}$.
A harmonic time-dependence of the fluctuating quantities is now assumed, such that 
$\hat {y} = \breve {y}{\text {e}}^{\text {i}\omega t}$, with 
$\omega$ the angular frequency, 
$\breve {y}$ the complex wave amplitude and 
$\text {i}^2={-1}$. Equation (2.9) is recast in the frequency domain as
\begin{equation} \dfrac{\text{d}}{\text{d}\kern0.06em x} \boldsymbol{\mathsf{I}}=\boldsymbol{\mathsf{A}}{(x, \:\omega, \bar{\dot{Q}})}\boldsymbol{\mathsf{I}} \quad \text{with} \ \boldsymbol{\mathsf{A}} ={-}{\left[\boldsymbol{\mathsf{E}}_{x}\right]}^{{-}1}\left(\dfrac{\text{i}\omega}{\bar{u}}\mathcal{I} + \dfrac{\left(\gamma-1\right)\bar{\dot{Q}}}{2\zeta \gamma \bar{p}\: \bar{u}} \boldsymbol{\mathsf{E}}_{s}\right), \end{equation}
where 
$\mathcal {I}$ is an identity matrix of order 3. Equation (2.11) is similar to the equation obtained by Duran & Moreau (Reference Duran and Moreau2013) but with an extra term, involving the matrix 
$\boldsymbol{\mathsf{E}}_{s}$, which accounts for mean heat transfer effects. The frequency and the heat source are specified in terms of the non-dimensional parameters 
$\varOmega = \omega L/\bar {c}_{t0}$ (axial Helmholtz number) and 
$\tilde {Q} = \bar {\dot {Q}}L/(\bar {p}_{t0}\; \bar {c}_{t0}),$ respectively, where 
$L$ is the length of the nozzle, and 
$\bar {p}_{t0}$ and 
$\bar {c}_{t0}$ correspond to the stagnation pressure and stagnation speed of sound at the inlet, respectively.
Equation (2.11) is solved using a Magnus-expansion-based method, that assumes the ansatz
\begin{equation} \boldsymbol{\mathsf{I}}(\xi, \varOmega, \tilde{Q})=[\exp [\boldsymbol{\mathsf{B}}(\xi, \varOmega, \tilde{Q})]] \boldsymbol{\mathsf{I}}_{0} \quad \text{with} \ \boldsymbol{\mathsf{B}}(\xi, \varOmega, \tilde{Q})=\sum_{k=0}^{\infty} \boldsymbol{\mathsf{B}}^{(k)}(\xi, \varOmega, \tilde{Q}), \end{equation}
with 
$\xi$ = 
$x/L$ denoting the dimensionless axial coordinate and 
$\boldsymbol{\mathsf{I}}_{0}$ the vector of invariants at the inlet. Here 
$\boldsymbol{\mathsf{B}}^{(k)}(\xi, \varOmega, \tilde {Q})$ represents the terms of the Magnus expansion each of order 
$O(\varOmega, \tilde {Q})^{k}.$ The terms in the expansion are obtained recursively as explained by Blanes et al. (Reference Blanes, Casas, Oteo and Ros2009). When the frequency of the fluctuating components (represented by 
$\varOmega$) and the heat exchange (represented by 
$\tilde {Q}$) are zero, the flow invariants defined in (2.8) are conserved and remain constant along the nozzle length. Any non-zero frequency and/or non-zero heat exchange results in a deviation which is captured by the higher order terms (
$k\neq 0$) in the Magnus expansion. The main novelty of this study is that, for the first time, a term capturing flow non-isentropicity is included as the expansion parameter along with frequency. To ensure convergence, the series may require the separate computation of transfer matrices for axial segments of the nozzle, which are subsequently multiplied for the final result (Blanes et al. Reference Blanes, Casas, Oteo and Ros2009). For example, for the isentropic subcritical mean nozzle flow considered in § 3, at low frequencies (
${\varOmega }/{(2{\rm \pi} )} \sim 0.0001\text {--}0.01$) the series exhibits fast convergence when applied to the entire nozzle length and does not require any axial segmentation. However, at higher frequencies, for example, when 
${{\varOmega }}/(2{\rm \pi} ) = 1$, the nozzle needs to be divided into at least 24 axial segments for fast convergence. This segmentation approach follows the fast convergence criterion given by
\begin{equation} \int_{0}^{\xi_{F}}\|\boldsymbol{A}(\xi)\|_{2} \,\mathrm{d} \xi<{\rm \pi}, \end{equation}
where 
$\xi _F$ determines the maximum length of the segments. The above (2.13) ensures faster convergence when satisfied, but may diverge if the value of the integral is larger than 
${\rm \pi}$ (Blanes et al. Reference Blanes, Casas, Oteo and Ros2009).
Finally, the flow invariants at the boundaries are transformed into three propagating waves (see figure 1): (i) a downstream-propagating acoustic wave, 
$w^+$, (ii) an upstream-propagating acoustic wave, 
$w^-,$ and (iii) an entropy wave, 
$w^s.$ These wave components at a particular location are represented using the wave vector 
$\boldsymbol{\mathsf{W}}=[w^+ \;w^-\; w^s]^{\rm T}$. In terms of the primitive variables, they read 
$w^+ = \breve {p} + \breve {u}M$, 
$w^- = \breve {p} - \breve {u}M,$ and 
$w^s = \breve {s}$, with 
$\breve {p}, \breve {u}$ and 
$\breve {s}$ denoting the complex Fourier coefficients for pressure, velocity and entropy, respectively.
The relation between the flow invariants and the wave components can then be obtained using (2.8) and is represented as 
$\boldsymbol{\mathsf{I}} = \boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{W}}$, where 
$\boldsymbol{\mathsf{D}}$ defines the relation between the wave components and the flow invariants. A transfer matrix 
$\boldsymbol{\mathsf{T}}$ is therefore defined to relate the wave vector at any location, 
$\xi$, in terms of the wave vector at the inlet (
$\boldsymbol{\mathsf{W}}_{0}$) as
\begin{equation} \boldsymbol{\mathsf{W}}_{\xi} = \boldsymbol{\mathsf{T}} \boldsymbol{\mathsf{W}}_{0}, \quad \text{where} \ \boldsymbol{\mathsf{T}} = \left[\boldsymbol{\mathsf{D}}_{\xi}\right]^{{-}1} \boldsymbol{\mathsf{C}} \; \boldsymbol{\mathsf{D}}_{0} \end{equation}
and 
$\boldsymbol{\mathsf{C}} = [\exp [\boldsymbol{\mathsf{B}}(\xi, \varOmega, \tilde {Q})]]$.
2.1. Subcritical nozzle flow configuration
For the case where the flow remains subsonic inside the nozzle, the acoustic system is determined by three external inputs, namely a downstream-propagating acoustic wave at the inlet, 
$w^+_{0, f}$, an upstream-propagating acoustic wave at the outlet, 
$w^-_{1, f},$ and an entropy wave at the inlet, 
$w^s_{0, f}$. Note that the subscripts ‘0’ and ‘1’ refer to the inlet and outlet of the nozzle, respectively, while ‘
$f$’ denotes that the wave is externally forced. We also have three waves as the outputs of the system, i.e. an upstream-propagating acoustic wave at the inlet, 
$w^-_0$, a downstream-propagating acoustic wave at the outlet, 
$w^+_{1}$ and an entropy wave at the outlet, 
$w^s_{1}$. An extended scattering matrix 
$\boldsymbol{\mathsf{S}}$ can then be defined to relate the incoming and outgoing wave vectors as
\begin{equation} \begin{bmatrix}
w^+_1 \\w^-_0 \\ w^s_1 \end{bmatrix} =
\boldsymbol{\mathsf{S}} \begin{bmatrix}
w^+_{0, f} \\w^-_{1, f} \\ w^s_{0, f} \end{bmatrix}.
\end{equation}
This matrix is obtained by algebraically rearranging the terms of the matrix 
$\boldsymbol{\mathsf{T}}$ in (2.14).
2.2. Supercritical nozzle flow without shocks
For a supercritical flow without any shock waves inside the nozzle, only two inputs, namely 
$w^+_{0, f}$ and 
$w^s_{0, f},$ can be applied as 
$w^-_{1}$ now corresponds to the slow downstream propagating wave at the nozzle outlet. The flow domain is then divided into subsonic (from the inlet 
$x_0$ to a location infinitesimally upstream of the throat 
$x_u$) and supersonic (from a location infinitesimally downstream of the throat 
$x_d$ to the outlet 
$x_1$) portions. Mass fluctuations cannot occur at the choked throat and this requirement is specified as a boundary condition (Marble & Candel Reference Marble and Candel1977) which takes the form 
${M}^{'}/M = 0$. At the location infinitesimally upstream of the throat, this gives
where 
$R_u$ and 
$R_s$ correspond to the acoustic and entropy reflection coefficients, respectively. The inlet acoustic and entropy forcing inputs, 
$w^+_{0, f}$ and 
$w^s_{0, f},$ respectively, along with (2.16) give the three input conditions. Equation (2.14) then takes the form
\begin{equation} \begin{bmatrix}
w^+_u \\ {w^-_{u, f}} \\ w^s_u
\end{bmatrix}\equiv\begin{bmatrix}
w^+_u \\ 0 \\ w^s_u
\end{bmatrix}=\underbrace{{\left[\boldsymbol{\mathsf{D}}_{1}^{'}\right]}^{{-}1}
\boldsymbol{\mathsf{C}}_{{sub}} \;
\boldsymbol{\mathsf{D}}_{0}}_{{\boldsymbol{\mathsf{T}}}^{'}}
\begin{bmatrix} w^+_{0, f}\\ w^-_{0} \\ w^s_{0,
f}\end{bmatrix}, \end{equation}
where matrix 
$\boldsymbol{\mathsf{D}}^{\prime}$ relates the wave components in terms of reflection coefficients to the flow invariants and matrix 
$\boldsymbol{\mathsf{C}}_{sub}$ connects these flow invariants at the inlet to the throat. Again rearranging the matrices with forcing inputs on one side and unknowns on the other gives,
\begin{equation} \begin{bmatrix}
w^+_u \\ w^-_0 \\ w^s_u \end{bmatrix} =
\boldsymbol{\mathsf{S}}_{s} \begin{bmatrix}
w^+_{0, f}\\ 0 \\ w^s_{0, f}\end{bmatrix} \Rightarrow
w^-_0 = S_s(2, 1)\;w^+_{0, f} + S_s(2, 3)\;w^s_{0,
f}, \end{equation}
where 
$\boldsymbol{\mathsf{S}}_{s}$ corresponds to the scattering matrix of the subsonic portion of the supercritical nozzle, given in terms of the elements of 
$\boldsymbol{\mathsf{T}}^{'}$. Separating the upstream-propagating wave at the inlet using (2.18) gives
\begin{equation} \begin{bmatrix}
w^+_u \\ 0 \\ w^s_u \end{bmatrix} =
{\boldsymbol{\mathsf{S}}'_{s}} \begin{bmatrix}
w^+_{0, f}\\ 0 \\ w^s_{0, f}\end{bmatrix} \quad
\text{where}\,\, {\boldsymbol{\mathsf{S}}'_{s}} =
\begin{bmatrix} S_s(1, 1)\quad 0 \quad S_s(1,
3)\\ 0 \quad 0 \quad 0 \\S_s(3, 1) \quad 0 \quad S_s(3, 3)
\end{bmatrix}. \end{equation}
The wave vector is unaltered at the throat and hence 
$\boldsymbol{\mathsf{W}}_{u} = \boldsymbol{\mathsf{W}}_d$. Note that 
$\boldsymbol{\mathsf{W}}_{u} = [w^+_{u} \; 0 \; w^s_{u}]^{\rm T}$ while 
$\boldsymbol{\mathsf{W}}_{d} = [w^+_{d} \; {w^-_{d}} \; w^s_{d}]^{\rm T}$ since 
$w^-_{u, f} = 0$. Denoting the transfer matrix in the supersonic portion of the nozzle by 
$\boldsymbol{\mathsf{T}}_{s}$, the wave vector relation 
$\boldsymbol{\mathsf{W}}_{1}= \boldsymbol{\mathsf{T}}_{s} \boldsymbol{\mathsf{W}}_{d}$ is obtained similar to (2.14), where 
$\boldsymbol{\mathsf{T}}_{s} = {[\boldsymbol{\mathsf{D}}_{1}]}^{-1} \boldsymbol{\mathsf{C}}_{{super}} \; \boldsymbol{\mathsf{D}}_{d}$ with matrix 
$\boldsymbol{\mathsf{D}}_d$ converting the wave components to flow invariants at the throat and matrix 
$\boldsymbol{\mathsf{C}}_{super}$ relating these invariants to the invariants at the outlet. This gives
\begin{equation} \begin{bmatrix}
w^+_1 \\ w^-_1 \\ w^s_1 \end{bmatrix}\; =
\boldsymbol{\mathsf{T}}_{s}
{\boldsymbol{\mathsf{S}}'_{s}} \begin{bmatrix}
w^+_{0, f}\\ 0 \\ w^s_{0, f}\end{bmatrix}.
\end{equation}
Equation (2.20) gives the fast (
$w^+_{1}$) and slow (
$w^-_{1}$) propagating acoustic waves and the entropy wave (
$w^s_{1}$) at the outlet as a function of acoustic (
$w^+_{0, f}$) and entropy (
$w^s_{0, f}$) forcing at the nozzle inlet. The upstream-propagating wave at the inlet (
$w^-_{0}$) in the subsonic portion is given by (2.18).
3. Results
In this section, we apply the presented theory to a converging–diverging nozzle (see figure 1) defined by the following area profile:
\begin{equation} \frac{A(x)}{A_{*}}=\left\{\begin{array}{@{}ll} \dfrac{1}{2}\left(\dfrac{A_{{0}}}{A_{*}}-1\right)\left[\cos \left({\rm \pi} \dfrac{x}{x_{*}}\right)+1\right]+1, & \text{if } x \in\left[0, x_{*}\right],\\ 1+\left(\dfrac{A_{{1}}}{A_{*}}-1\right) \dfrac{x-x_{*}}{L-x_{*}}, & \text{if } x \in\left[x_{*}, L\right], \end{array}\right. \end{equation}
where 
$A_* = 0.002\,\text {m}^2$ corresponds to the throat area (at the throat location 
$x_* = 0.15$ m), 
$A_{0}/{A_*} = 2.1$ and 
$A_{1}/{A_*} = 1.18$. The gas constant and adiabatic index are taken to be 
${R}_{g}=287\,\text {J}\,\text {kg}^{-1}\,\text {K}^{-1}$ and 
$\gamma =1.4,$ respectively.
Two different flow cases are considered: (i) subcritical and (ii) supercritical (without shocks). For both flow cases, numerical solutions of the mean flow are obtained for different levels of heat transfer, set by 
$\tilde {Q}.$ In the following, we assume for simplicity 
$\tilde {Q}$ to be constant axially, but the approach can assume any spatial distribution. The subcritical case is defined by an inlet Mach number of 
$M_{0} = 0.2.$ The supercritical case assumes 
$M_{0} = 0.29$ for the isentropic case. When heat is added, the inlet Mach number ranges from 
$M_0=0.28$ to 
$M_0=0.30$ for 
$\tilde {Q}=0.3$ and 
$\tilde {Q}=-0.5,$ respectively. The mean flow for the subcritical case is obtained by solving the first-order system of differential equations of (2.3) using higher order implicit schemes: for example, a fourth-order Runge–Kutta was used here. For the supercritical case, a modified finite-difference-based MacCormack scheme is used to compute the mean flow. The evolution of the Mach number and mean temperature inside the nozzle are presented in figure 2. As expected, the temperature at the outlet is higher when heat is added than for the isentropic flow. For the Mach number distribution, we observe different trends for the subcritical and supercritical cases. For the subcritical nozzle flow, the Mach number at the outlet increases when heat is added owing to the decrease of density and subsequent increase of the flow velocity to satisfy continuity. For the supercritical nozzle flow, adding heat reduces the Mach number at the outlet, explained by the increase of the speed of sound due to the temperature rise.
Mach number, 
$M$, (lines) and mean temperature, 
$\bar {T}$, (lines with markers) for (
$a$) subcritical case with 
$\tilde {Q}$ values of 
$-0.5$ (thick dash-dot-dash line), 0 (—–), 0.5 (dashed line); (
$b$) supercritical case for 
$\tilde {Q}$ values of 
$-0.5$ (thick dash-dot-dash line), 0 (—–), 0.3 (dashed line).
3.1. Effect of heat transfer at zero frequency
In this section, we explore the influence of heat transfer on the unsteady response of the nozzle at zero frequency (
$\varOmega = 0)$. The heat source/sink is varied in the range 
$\tilde {Q}\in [-0.5, 0.5]$ for the subcritical flow and 
$\tilde {Q}\in [-0.5, 0.3]$ for the supercritical flow. For 
$\tilde {Q}\gtrsim 0.3$ in the supercritical regime, a normal shock appears. While the model can be extended to nozzles sustaining normal shocks (Duran & Moreau Reference Duran and Moreau2013), the current validation is limited to shock-free nozzles only. Figures 3 and 4 show the coefficients of the acoustic transfer functions for the subcritical and supercritical nozzle flow cases, respectively. The order of the Magnus expansion considered is 
$k=5$. The model predictions are compared to numerical solutions and an excellent agreement is observed.
Transfer functions for the subcritical nozzle flow at zero frequency predicted by the present model (∗, blue), the compact isentropic model of Marble & Candel (Reference Marble and Candel1977) (—, red) and numerical solutions (dotted line).
Transfer functions for the supercritical nozzle flow at zero frequency predicted by the present model (∗, blue), the compact isentropic model of Marble & Candel (Reference Marble and Candel1977) (—, red) and numerical solutions (dotted line).
To assess the importance of the heat transfer in the acoustic solution, the isentropic and compact theory of Marble & Candel (Reference Marble and Candel1977) is also plotted. This model requires only two parameters, namely the Mach number at the inlet and at the outlet. Here, we feed those parameters to the Marble and Candel model based on the mean flow numerical solution with heat transfer. This results in predictions where the perturbations are isentropic, but the mean flow is corrected to account for non-isentropicity. As observed, the effect of heat transfer on the results is strong: the trends for some of the coefficients are reversed when compared with isentropic theory and the relative errors are as high as several hundred percent. The differences can be explained based on the governing equations of the fluctuating variables. The isentropic and non-isentropic equations only differ in the source term owing to heat transfer on the right-hand side of (2.4) and (2.6). These terms represent sources of fluctuating mass flow rate, stagnation temperature and entropy that are produced by the interaction of acoustic waves with the steady heat transfer. The variations of fluctuating mass flow rate and stagnation temperature act as an additional source of sound. The variations of entropy, on the other hand, act as a dipole source of sound through the momentum equation (2.5). For the subcritical case, the dominant effect is the former. This can be argued based on the relatively low values of the acoustic response to entropy fluctuations (figure 3b,e). For values of around 
$\tilde {Q}\approx -0.4,$ this is especially notable, as the sound produced by entropy tends to zero, but the acoustic coefficients substantially differ from the predictions of isentropic theory. For the supercritical configuration, on the other hand, both effects are similarly important, as can be deduced from the high values of entropy-sound generation (figure 4a,b). An interesting trend is observed for the slow-propagating downstream wave, 
$w^-_1,$ for the supercritical case: for values 
$0.1\lesssim \tilde {Q} \lesssim 0.3,$ the coefficient exhibits a strong sensitivity when forced by either upstream acoustic waves or entropy waves. This effect, not observed for the fast-propagating downstream wave, may be linked to the fact that the Mach number is close to unity along the whole diverging portion of the nozzle in this range and, therefore, the propagation speed of that wave is close to zero.
3.2. Effect of heat transfer for non-zero frequencies
We now turn our attention to the effect of heat transfer at non-zero frequencies. Figure 5 shows the variation of magnitude and phase of the acoustic transfer functions for the subcritical nozzle flow, for both isentropic (
$\tilde {Q}=0$) and non-isentropic (
$\tilde {Q}=-0.5$) cases. It can again be observed that the model estimates closely match the numerical predictions. Similar trends are observed for all of the coefficients in the frequency domain. On comparing the isentropic and non-isentropic results for the magnitude of some of the coefficients, i.e. 
$w^+_1/w^+_{0,f},$ 
$w^+_1/w^-_{1,f}$ and 
$w^-_0/w^-_{1,f}$, an offset was observed at zero frequency which is carried to higher frequencies. For the rest of the coefficients, large differences exist at low frequencies, but quickly disappear with increasing frequencies. The phases for all of the coefficients are remarkably similar, with some differences observed for the coefficients defined in the downstream duct. This difference is expected due to the large difference in the flow Mach numbers and speed of sound in the divergent portion of the nozzle for nozzles with different levels of 
$\tilde {Q}$. This leads to different propagation velocities of the acoustic waves in the diverging portion of the nozzle and therefore to the differences in phase.
Transfer functions for the subcritical nozzle flow as a function of the frequency 
$\varOmega$. Numerical (dotted line) and model solutions (
$\bigcirc$) for 
$\tilde {Q} = 0$. Numerical (blue dashed line) and model solutions (
${\lozenge }$, blue) for 
$\tilde {Q} = -0.5$. Marble & Candel (Reference Marble and Candel1977) compact solution 
$(\bigstar )$.
Figure 6 shows the variation of the acoustic transfer functions as a function of the frequency for the supercritical nozzle flow with 
$\tilde {Q} = -0.5$ and 
$0$. The model estimates are accurate for both cases. The trends obtained for the reflection coefficient in the upstream duct 
$w^-_0$ are very similar to those obtained for the subcritical case and have been omitted for the sake of brevity. This is consistent with 
$w^-_0$ not depending on the supersonic portion of the nozzle such that it can be obtained using (2.15) alone. The acoustic reflection at the outlet (
$w^-_{1}$) has a low-pass-like behaviour for both entropy (
$w^s_{0, f}$) and acoustic (
$w^+_{0, f}$) forcing at the inlet.
Transfer functions for the supercritical nozzle flow as a function of the frequency 
$\varOmega$. Numerical (dotted line) and model solutions (
$\bigcirc$) for 
$\tilde {Q} = 0$. Numerical (blue dashed line) and model solutions (
${\lozenge }$, blue) for 
$\tilde {Q} = -0.5$. Marble & Candel (Reference Marble and Candel1977) compact solution 
$(\bigstar )$.
4. Conclusions
This work proposed a model for the acoustic and entropic transfer functions of non-isentropic nozzle flows where the non-isentropicity arises due to a steady source/sink of heat. Physically, this source term can be produced by heat exchange with surroundings or by chemical reactions. The flow was modelled as a quasi-one-dimensional and inviscid flow. A solution method was proposed based on the Magnus expansion, which does not impose any limitations on the frequency or Mach number that can be considered. In this work, a measure of the flow non-isentropicity has been used as an expansion parameter for the first time. The solution was successfully validated against numerical simulations of the quasi-one-dimensional linearised Euler equations in a converging–diverging nozzle for different heat transfer and frequency levels. Two types of nozzle flow were considered: subcritical and supercritical (without shocks). For the acoustically compact case, it was shown that applying the Marble & Candel (Reference Marble and Candel1977) solution in the presence of heat transfer can lead to large prediction errors, even when the corrected mean flow Mach numbers are provided. The heat transfer was found to strongly affect the transfer functions of the nozzle, with some acoustic transfer functions being amplified and others attenuated, for increasing heat transfer.
Funding
This work was supported by the European Research Council (ERC) Consolidator Grant AFIRMATIVE (2018-2023) and the Inlaks Shivdasani Foundation.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Linearised Euler equations in terms of flow invariants
Using the definition of flow invariants from (2.8), we can write
\begin{equation} \overbrace{\left[\begin{array}{c} I_A \\ I_B \\ I_C \end{array}\right]}^{\boldsymbol{\mathsf{I}}} =\overbrace{\left[\begin{array}{cccc} 1 & 1 & -1 \\ \dfrac{\gamma-1}{\zeta} & \dfrac{(\gamma-1)M^2}{\zeta} & \dfrac{1}{\zeta} \\ 0 & 0 & 1 \\ \end{array}\right]}^{\boldsymbol{\mathsf{D}}_{I}^{P}} \overbrace{\left[\begin{array}{c} \hat{p} \\ \hat{u} \\ \hat{s} \end{array}\right]}^{\boldsymbol{\mathsf{P}}}, \quad \Rightarrow \boldsymbol{\mathsf{P}}=\left[\boldsymbol{\mathsf{D}}_{I}^{P}\right]^{{-}1} I. \end{equation}By using the time-averaged mean flow conservation equations (2.3), we can write
\begin{equation} \dfrac{\text{d}\bar{u}}{\text{d}\kern0.06em x} = \dfrac{\bar{u}}{\zeta}\left(\dfrac{1}{M}\dfrac{\text{d}M}{\text{d}\kern0.06em x} + \dfrac{\left(\gamma-1\right)\bar{\dot{Q}}}{2 \gamma \bar{p}\:\bar{u}}\right). \end{equation}(i) Invariant formulation for normalised fluctuating mass,
$I_A$:Adding and subtracting (2.5) and (2.6), respectively, from (2.4) and using (A2) gives,
(A 3)Using (2.5) we substitute for${ \partial \hat {u}}/{ \partial t }$ and (A2) is used to replace ${\text {d}\bar {u}}/{\text {d}{x}}$. This gives
(A 4)where\begin{equation} \dfrac{ \text{D} I_A}{ \text{D} t } + \dfrac{\bar{u}}{\zeta} \left(\dfrac{1}{M}\dfrac{\text{d}M}{\text{d}\kern0.06em x} + \dfrac{\left(\gamma-1\right)\bar{\dot{Q}}}{2 \gamma \bar{p}\:\bar{u}}\right)\left(2\hat{u}-\left(\gamma-1\right) \hat{p}-\hat{s}\right) + \bar{u}\dfrac{\partial \hat{u}}{\partial x} + \dfrac{\bar{c}^2}{\bar{u}}\dfrac{\partial \hat{p}}{\partial x} = 0, \end{equation}(A 5)Taking the factor\begin{equation} \dfrac{{\text{D}} {(\;)}}{{\text{D}} t} = \dfrac{\partial(\;)}{\partial t} + \bar{u}\dfrac{\partial(\;)}{\partial x}. \end{equation}$1/M^2$ and writing out all the terms explicitly results in
(A 6)Algebraic rearrangement simplifies the above equation to$$\begin{gather} \dfrac{ \text{D} I_A}{ \text{D} t } + \dfrac{1}{M^2}\left(2M\dfrac{\text{d}M}{\text{d}x}\dfrac{\bar{u}}{\zeta} \hat{u} - \left(\gamma - 1\right)M\dfrac{\text{d}M}{\text{d}x}\dfrac{\bar{u}}{\zeta}\hat{p} - M\dfrac{\text{d}M}{\text{d}x} \dfrac{\bar{u}}{\zeta}\hat{s}\right.\nonumber\\ \left.+ \dfrac{\left(\gamma-1\right)M^2\bar{\dot{Q}}}{2 \gamma \bar{p}\:\bar{u}}\left(2\hat{u}-\left(\gamma-1\right) \hat{p}-\hat{s}\right) + \bar{u}\dfrac{\partial \left(M^2 \hat{u}\right)}{\partial x} - 2M\dfrac{\text{d}M}{\text{d}x}\bar{u} \hat{u} + \bar{u}\dfrac{\partial \hat{p}}{\partial x} \right) = 0. \end{gather}$$(A 7)Using (2.8) we replace$$\begin{gather} \dfrac{ \text{D} I_A}{ \text{D} t } + \dfrac{1}{M^2}\left( - \left(\gamma - 1\right)M\dfrac{\text{d}M}{\text{d}x}\dfrac{\bar{u}}{\zeta}\left( \hat{u}M^2 + \hat{p} + \dfrac{\hat{s}}{\gamma-1}\right) + \bar{u}\dfrac{\partial}{\partial x}\left( \hat{u}M^2 + \hat{p} + \dfrac{\hat{s}}{\gamma-1}\right)\right. \nonumber\\ \left.+\dfrac{\bar{u}}{\gamma - 1} \dfrac{\partial \hat{s}}{\partial x} + \dfrac{\left(\gamma-1\right)M^2\bar{\dot{Q}}}{2 \gamma \bar{p}\:\bar{u}}\left(2\hat{u}-\left(\gamma-1\right) \hat{p}-\hat{s}\right) \right) = 0. \end{gather}$$$( \hat {u}M^2 + \hat {p} +{\hat {s}}/({\gamma -1}))$ with ${\zeta I_B}/({\gamma - 1})$ and $\hat {s}$ with $I_C$. This gives
(A 8)On further carrying out the differential$$\begin{gather} \dfrac{ \text{D} I_A}{ \text{D} t } + \dfrac{1}{M^2}\left( - M\dfrac{\text{d}M}{\text{d}x}\bar{u} I_B + \bar{u}\dfrac{\partial}{\partial x}\left( \dfrac{\zeta I_B}{\gamma - 1}\right)+ \dfrac{\bar{u}}{\gamma - 1} \dfrac{\partial I_C}{\partial x}\right) \nonumber\\ + \dfrac{\left(\gamma-1\right)\bar{\dot{Q}}}{2 \gamma \bar{p}\:\bar{u}}\left(2\hat{u}-\left(\gamma-1\right) \hat{p}-\hat{s}\right) = 0. \end{gather}$$${\partial }/{\partial x}({\zeta I_B}/({\gamma - 1}))$ in the above equation, we obtain
(A 9)\begin{equation} \dfrac{ \text{D} I_A}{ \text{D} t } + \dfrac{1}{M^2}\dfrac{\bar{u}}{\gamma-1}\left(\zeta\dfrac{\partial I_B}{\partial x} - \dfrac{\partial I_C}{\partial x} \right) + \dfrac{\left(\gamma-1\right)\bar{\dot{Q}}}{2 \zeta \gamma \bar{p}}\left(2\hat{u}-\left(\gamma-1\right) \hat{p}-\hat{s}\right) = 0. \end{equation}(ii) Invariant formulation for normalised stagnation temperature fluctuations,
$I_B$:On multiplying (2.5) by
$M^2$, and using (A2), we can write
(A 10)Again using (2.8), we can replace$$\begin{gather} \dfrac{ \partial}{ \partial t }\left(\hat{u}M^2 \!+ \hat{p} + \dfrac{\hat{s}}{\gamma\!-\!1}\right) + \bar{u}\dfrac{ \partial}{ \partial x }\left(\hat{u}M^2 + \hat{p} + \dfrac{\hat{s}}{\gamma\!-\!1}\right) -\dfrac{1}{\gamma\!-\!1}\dfrac{\text{D}\hat{s}}{\text{D} t} - \dfrac{ \partial \hat{p}}{ \partial t } -2\bar{u}M\dfrac{\text{d}M}{\text{d}x}\hat{u}\nonumber\\ +\dfrac{\bar{u}}{\zeta}\left(M\dfrac{\text{d}M}{\text{d}\kern0.06em x} + \dfrac{\left(\gamma-1\right)\bar{\dot{Q}}}{2 \gamma \bar{p}\:\bar{u}}\right)\left(2\hat{u}-\left(\gamma-1\right) \hat{p}-\hat{s}\right) = 0. \end{gather}$$$( \hat {u}M^2 + \hat {p} +{\hat {s}}/({\gamma -1}))$ by ${\zeta I_B}/({\gamma - 1})$. Here ${\text {D}\hat {s}}/{\text {D}t}$ and ${\partial \hat {p}}/{\partial t}$ are also substituted using (2.6) and (2.4), respectively. This gives
(A 11)Further using (2.8), we replace$\hat {p} +\hat {u}$ with $I_A + I_C$. This results in
(A 12)On multiplying the resultant (A12) with$(\gamma -1)/\zeta$ and rearranging, we get
(A 13)$$\begin{gather} \dfrac{ \text{D} I_B}{ \text{D} t } + \dfrac{\left(\gamma-1\right)}{\zeta} {\bar{u}} \dfrac{\partial \left(I_{A} +I_C\right)}{\partial x} + \dfrac{\left(\gamma-1\right)\bar{\dot{Q}}}{\gamma \bar{p} \zeta} \left(\hat{p}\left(\gamma^2-\dfrac{{\left(\gamma-1\right)}^{2} M^2 }{2\zeta}\right)\right.\nonumber\\ \left.+ \hat{u} \left(\gamma + \dfrac{\left(\gamma-1\right)M^2}{\zeta}\right) - \hat{s} \dfrac{\left(\gamma-1\right)M^2}{2\zeta}\right)= 0. \end{gather}$$(iii) Invariant formulation for normalised fluctuating entropy,
$I_C$:Using
$I_C = \hat {s}$, we also have from (2.6)
(A 14)\begin{equation} \dfrac{ \text{D} I_C}{ \text{D} t } + \dfrac{\left(\gamma-1\right)\bar{\dot{Q}}}{\gamma \bar{p}}\left(\gamma\hat{p} + \hat{u}\right) = 0. \end{equation}
The above (A9), (A13) and (A14) can be written in a simplified matrix form as,
\begin{equation} \dfrac{\partial}{\partial t} \boldsymbol{\mathsf{I}} +\bar{u} \boldsymbol{\mathsf{E}}_{x} \frac{\partial}{\partial x} \boldsymbol{\mathsf{I}}+ \dfrac{\left(\gamma -1\right) \bar{\dot{Q}}}{2\zeta \gamma \bar{p}}\boldsymbol{\mathsf{E}}_{s, primitive} \boldsymbol{\mathsf{P}} = 0, \end{equation}where
\begin{equation} \boldsymbol{\mathsf{E}}_{s, primitive}=\left[\begin{array}{ccc} {-\left(\gamma-1\right)} & 2 & -1 \\ 2\gamma^2-\dfrac{{\left(\gamma-1\right)}^{2} M^2 }{\zeta} & 2\gamma + \dfrac{2\left(\gamma-1\right)M^2}{\zeta} & -\dfrac{\left(\gamma-1\right)M^2}{\zeta}\\ 2\zeta \gamma & 2\zeta & 0 \end{array}\right], \end{equation}
and 
$\boldsymbol{\mathsf{E}}_x$ is given in (2.9).
Using (A1), we can write
\begin{equation} \boldsymbol{\mathsf{E}}_{s, primitive} \boldsymbol{\mathsf{P}} \equiv \boldsymbol{\mathsf{E}}_{s, primitive} \left[\boldsymbol{D}_{I}^{P}\right]^{{-}1} \boldsymbol{\mathsf{I}} = \boldsymbol{\mathsf{E}}_s \boldsymbol{\mathsf{I}}, \end{equation}
where 
$\boldsymbol{\mathsf{E}}_s$ is as given by (2.10). Thus the final system of linearised Euler equations in terms of flow invariants takes the form as in (2.9).
Open access
