1. Introduction
Brass instruments are understood to sound brassy due to nonlinear wave steepening within the instrument (Hirschberg et al. Reference Hirschberg, Gilbert, Msallam and Wijnands1996). This has been observed experimentally for a trombone (Hirschberg et al. Reference Hirschberg, Gilbert, Msallam and Wijnands1996; Rendón et al. Reference Rendón, Orduña Bustamante, Narezo, Pérez-López and Sorrentini2010) and for a trumpet (Pandya, Settles & Miller Reference Pandya, Settles and Miller2003; Rendón et al. Reference Rendón, Ezeta and Pérez-López2013), and is supported by relatively simple physical models (e.g. Gilbert, Menguy & Campbell Reference Gilbert, Menguy and Campbell2008). In these studies, attention is often given to the width, or bore, of an instrument, and its variation along its length. Varying the bore so that it progressively gets wider lowers the amplitude of the outward-propagating wave, resulting in less nonlinear steepening and a less brassy sound. Conversely, a cylindrical bore that does not widen (such as is necessary for a trombone slide to be able to slide) allows nonlinear steepening to occur, exciting higher harmonics and leading to a more brassy sound. Myers et al. (Reference Myers, Pyle, Gilbert, Campbell, Chick and Logie2012) even proposed a dimensionless brassiness parameter
$B$
based on the bore variation in order to quantitatively classify musical instruments by their brassiness. Most of the modelling in this area assumes plane-wave propagation and straight ducts, although Fernando & Druon (Reference Fernando and Druon2011) allowed for non-plane-wave propagation by introducing a cross-duct modal expansion while retaining the straight-duct approximation.
However, very few brassy-sounding musical instruments are straight; for example, a traditional B
$\flat$
trumpet has a duct that, if straightened, would be approximately three times as long as the instrument. In order to investigate the acoustics of curved ducts, Félix & Pagneux (Reference Félix and Pagneux2001) generalised the multi-modal method previously developed by Pagneux, Amir & Kergomard (Reference Pagneux, Amir and Kergomard1996) for ducts with width variation. While initially developed in two dimensions, it was subsequently extended to three (Félix & Pagneux Reference Félix and Pagneux2002). The multi-modal method involves the projection of the curved-duct acoustics equations for pressure and velocity onto a basis of straight-duct modes, converting the governing partial differential equations (PDEs) into an infinite coupled set of ordinary differential equations (ODEs) for the amplitude of each mode. The ratio of the velocity and pressure, known as the admittance, is then taken, and a Riccati-style equation for the admittance is solved. As well as being computationally tractable, this formulation has the potential to provide better physical intuition than, for example, direct numerical simulation, as the admittance encodes the downstream effects of the duct geometry independently of the form of the acoustic wave introduced upstream. This will be seen below, as we will follow a similar procedure here. However, this multi-modal method depends on the linearity of the acoustic equations, and so is not able to model the weak nonlinearity necessary for investigating brassiness.
The tuning of a wind instrument is given by the effective length of the duct, and acoustic modes travelling around the outside of a duct bend would see a different effective length compared with acoustic modes travelling around the inside of a duct bend, resulting in a different resonant pitch. Brass instruments shift their pitch slightly when played louder. Could perhaps nonlinear steepening affect the effective length of a duct bend? Could instruments be designed for pitch stability as the sound amplitude is varied? Or, if curved ducts behave similarly to straight ducts, is it possible to say why and to quantify the extent to which they do so? Such an analysis would require a framework capable of analysing nonlinear acoustics in curved ducts.
Recently, McTavish & Brambley (Reference McTavish and Brambley2019) combined the curved-duct linear multi-modal method of Félix & Pagneux (Reference Félix and Pagneux2001) with the straight-duct nonlinear multimodal method of Fernando & Druon (Reference Fernando and Druon2011) to describe the weakly nonlinear acoustics of curved ducts in two dimensions. Even restricted to two dimensions, the algebraic complexity of the resulting system was prohibitive, and so a linear matrix and nonlinear tensor convolution notation was developed. Because the matrices and tensors involved in the governing equations varied as the duct geometry varied, the technique was computationally inefficient; nonetheless, it was possible to analyse a range of curved geometries in both the linear and nonlinear regimes.
The aim of the work presented here was originally to extend the two-dimensional analysis of McTavish & Brambley (Reference McTavish and Brambley2019) to three dimensions. In so doing, a modified mathematical framework was developed that combines both two- and three-dimensional cases, is computationally more efficient and is easier to understand in terms of the interplay between the nonlinearity and the duct geometry (consisting of width variation, curvature and, in three dimensions, torsion). Section 2 derives this mathematical framework from the governing fluid mechanics equations, resulting in the multimodal set of governing ODEs valid in both two and three dimensions given in § 2.6. These equations are then solved using the multi-modal method in § 3, first by introducing the admittance (§ 3.1) and then by truncating and numerically solving (§ 3.5). A number of numerical results are then presented in § 4, including validation for a straight duct (§ 4.1), a constant curvature bend in two and three dimensions (§§ 4.2 and 4.3) and an exponential horn (§ 4.4), along with more interesting cases involving combinations of curvature, width variation, torsion and nonlinearity. In particular, the leakage of acoustics from a duct owing to either curvature or nonlinearity is demonstrated in § 4.6, a direct comparison of two-dimensional results against three-dimensional results made possible by the combined mathematical framework is given in § 4.8 and an investigation of the effective acoustic length of a bend is given in § 4.9. The results are summarised and discussed in the concluding § 5, together with possibilities for future work. Matlab source code that produces the results shown here is provided in the supplementary material available at https://doi.org/10.1017/jfm.2026.11220, along with videos animating some of the figures.
2. Formulation of the governing equations
In this section, we derive the governing equations in the weakly nonlinear limit that will be solved in subsequent sections. The derivation follows that of McTavish & Brambley (Reference McTavish and Brambley2019), although here we do not limit ourselves to only two dimensions.
2.1. Derivation of the weakly nonlinear equations
We start with the equations of mass and momentum conservation for an inviscid fluid
\begin{align} \frac {\partial \hat {\rho }}{\partial \hat {t}} + \boldsymbol{\hat{\nabla }\boldsymbol{\cdot }}(\hat {\rho }\boldsymbol{\hat {u}}) = 0, \quad \hat {\rho }\! \left (\frac {\partial \hat {\boldsymbol{u}}}{\partial \hat {t}} + \boldsymbol{\hat {u}\boldsymbol{\cdot }\hat {\boldsymbol{\nabla }}\hat {u}}\right ) = -\hat {\boldsymbol{\nabla }}\!\hat {p},\end{align}
where
$\hat{\rho}$
,
$\hat{p}$
and
$\hat{\boldsymbol{u}}$
are the fluid density, pressure, and velocity respectively. Dimensional quantities are written with hats here. In this paper, we work in either two or three dimensions; note that in
$n$
dimensions, these give us
$n+1$
equations in
$n+2$
variables (one density,
$n$
velocities and one pressure). We therefore close the system of equations using the equation of state in combination with the thermal energy equation for an adiabatic gas (consistent with an inviscid non-heat-conducting fluid)
\begin{align} \hat {p}=\hat {p}(\hat {S},\hat {\rho }), \quad \frac {{\rm{D}}\hat {S}}{{\rm{D}}\hat {t}} = 0,\end{align}
where
$\hat{S}$
is the fluid entropy. Since we are interested in acoustics, we expand these equations for a small perturbation about an ambient state of rest. To do so, we first need some notion of size, so we introduce reference quantities and non-dimensionalise. Let the ambient density, pressure and entropy be denoted by
$\hat {\rho }_0$
,
$\hat {p}_0$
and
$\hat {S}_0$
, respectively. The ambient speed of sound
$\hat {c}_0$
is then given by
\begin{equation} \hat {c}_0^2 = \frac {\partial \hat {p}}{\partial \hat {\rho }}\bigg |_{\hat {S}}(\hat {\rho }_0,\hat {S}_0), \end{equation}
where the subscript zero here and elsewhere in this section denotes evaluation at the ambient state. Using
$\hat {\rho }_0$
and
$\hat {c}_0$
, non-dimensional variables are defined by setting
\begin{align} &\hat {\rho } = \hat {\rho }_0\rho , \quad \hat {\boldsymbol{u}} = \hat {c}_0\boldsymbol{u}, \quad \hat {p} = \hat {\rho }_0 \hat {c}_0^2 p. \end{align}
We also need spatial and temporal reference scales. The spatial reference length scale is denoted
$\hat {\ell }_0$
, typically given by the inlet radius of the duct. The time scale is then
$\hat {c}_0/\hat {\ell }_0$
, so the non-dimensional operators are
\begin{align} &\hat {\boldsymbol{\nabla }} = \frac {1}{\hat {\ell }_0}\boldsymbol{\nabla } \quad {\rm and} \quad \frac {\partial }{\partial \hat {t}} = \frac {\hat {c}_0}{\hat {\ell }_0}\frac {\partial }{\partial t}. \end{align}
We now perturb about an ambient state of rest
\begin{align} &\rho = 1 + \rho^{\prime}, \quad \boldsymbol{u} = \boldsymbol{0} + \boldsymbol{u}', \quad p = p_0 + p', \end{align}
taking
$\rho^{\prime} \sim p' \sim |\boldsymbol{u}'| \sim M$
, where
$M\ll 1$
is the perturbation Mach number. Since we are interested in weak nonlinearity, we will neglect terms of
$O(M^3)$
or smaller but retain both the linear acoustic
$O(M)$
and weakly nonlinear
$O(M^2)$
terms. Using the non-dimensionalisations (2.4, 2.5) in the governing (2.1) and substituting the perturbed quantities (2.6) then results in
\begin{align} \frac {\partial \rho^{\prime}}{\partial t} + \boldsymbol{\nabla \boldsymbol{\cdot }u}' &= -\rho^{\prime}\boldsymbol{\nabla \boldsymbol{\cdot }u}' - \boldsymbol{u}'\boldsymbol{\cdot \boldsymbol{\nabla }}\rho^{\prime} = \rho^{\prime}\frac {\partial \rho^{\prime}}{\partial t} - \boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla \!}\rho^{\prime} + O(M^3), \end{align}
\begin{align} \frac {\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{\nabla \!} p' &= -\boldsymbol{u}'\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}' - \rho^{\prime}\frac {\partial \boldsymbol{u}'}{\partial t} + O(M^3). \end{align}
Here, and in what follows, we write the linear acoustic
$O(M)$
quantities on the left-hand side of equations and the nonlinear
$O(M^2)$
terms on the right-hand side. Note that in (2.7a
) the linear expression
$\boldsymbol{\nabla \boldsymbol{\cdot }u}' = -\partial \rho^{\prime}/\partial t + O(M^2)$
was substituted in the right-hand side, resulting in an expression that is still correct to
$O(M^3)$
. This is a technique we will make further use of below.
Since we are primarily interested in acoustic perturbations, we will neglect entropy perturbations, so that
$S' \equiv 0$
and so
$\hat {S} \equiv \hat {S}_0$
everywhere. This means all perturbations are adiabatic, and so we may eliminate either the pressure or density perturbation using the expanded equation of state
\begin{align} p' = \rho^{\prime} + \frac {1}{2}\frac {\partial ^2 p}{\partial \rho ^2}\bigg |_{S}{\rho ^{\prime 2}} + O(M^3) \quad \Leftrightarrow \quad \rho^{\prime} = p' - \frac {1}{2}\frac {\partial ^2 p}{\partial \rho ^2}\bigg |_{S}{p}^{\prime 2} + O(M^3), \end{align}
where the second-order partial derivative is evaluated at the ambient state
$(\rho _0, S_0)$
, and the identity
$p' = \rho^{\prime} + O(M^2)$
has been used to rearrange between the left and right expressions. If we consider a perfect gas, a consequence of adiabaticity is that
$(\hat {p}/\hat {p}_0)=(\hat {\rho }/\hat {\rho }_0)^{ \unicode{x03B3}}$
, where the adiabatic index
$ \unicode{x03B3}$
is the ratio of specific heats. Consequently, for a perfect gas
\begin{align} \hat {c}_0^2 = \frac {\partial \hat {p}}{\partial \hat {\rho }}\bigg |_{\hat {S}} = \frac { \unicode{x03B3}\hat {p}_0}{\hat {\rho }_0}, \quad {\rm and} \quad \frac {\partial ^2 \hat {p}}{\partial \hat {\rho }^2}\bigg |_{\hat {S}} = \unicode{x03B3}( \unicode{x03B3}-1)\frac {\hat {p}_0}{\hat {\rho }_0^2} = ( \unicode{x03B3}-1)\frac {\hat {c}_0^2}{\hat {\rho }_0}, \end{align}
so that
$p_0 = 1/ \unicode{x03B3}$
and
$\partial ^2 p/\partial \rho ^2|_S = ( \unicode{x03B3}-1)$
when
$\rho _0=1$
and
$c_0=1$
. In what follows, for simplicity we will use this perfect gas notation, although the derivation is general provided it is taken that
$ \unicode{x03B3} = 1+\partial ^2 p/\partial \rho ^2|_S$
and
$p_0$
is not assumed to be
$1/ \unicode{x03B3}$
.
We now eliminate the density perturbation
$\rho^{\prime}$
using (2.8), and, correct to
$O(M^2)$
and dropping the
$O(M^3)$
from each equation for brevity, this finally results in
\begin{align} \frac {\partial p'}{\partial t} + \boldsymbol{\nabla \boldsymbol{\cdot }u}' &= \unicode{x03B3} p'\frac {\partial p'}{\partial t} - \boldsymbol{u}'\boldsymbol{\cdot \boldsymbol{\nabla }}p' = \frac {1}{2}\frac {\partial }{\partial t}\big ( \unicode{x03B3} {p}^{\prime 2} + |\boldsymbol{u}'|^2\big ),\\[-12pt]\nonumber \end{align}
\begin{align} \frac {\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{\nabla \!} p' &= -\boldsymbol{u}'\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}' - p'\frac {\partial \boldsymbol{u}'}{\partial t} = \frac {1}{2}\boldsymbol{\nabla }({p}^{\prime 2}) - \boldsymbol{u}'\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}', \end{align}
where again the
$O(M)$
identity
$\partial \boldsymbol{u'}/\partial t = -\boldsymbol{\nabla }p' + O(M^2)$
from the left-hand side of (2.10b
) has been used to rearrange the right-hand sides of both equations.
2.2. Further assumptions
Equations (2.10) may be simplified further. Taking
$\boldsymbol{\nabla }$
of both sides of (2.10b
), we apply symmetry of mixed partial derivatives on
$\boldsymbol{\nabla }\boldsymbol{\nabla \!} p'$
and find that the spin tensor of the velocity is linearly constant in time, i.e.
\begin{equation} \frac {\partial }{\partial t}\big (\boldsymbol{\nabla }\boldsymbol{u}'- (\boldsymbol{\nabla }\boldsymbol{u}')^{{T}}\big ) = 0 + O(M^2). \end{equation}
If we assume all acoustic variables to be periodic, we can define a time-average
$\langle \bullet \rangle$
, and then deduce that
\begin{equation} \boldsymbol{\nabla }\boldsymbol{u}' - (\boldsymbol{\nabla }\boldsymbol{u}')^{{T}} = \big \langle \boldsymbol{\nabla }\boldsymbol{u}' - (\boldsymbol{\nabla }\boldsymbol{u}')^{{T}}\big \rangle + O(M^2). \end{equation}
Having worked to a relatively high degree of generality so far, we now outline a hierarchy of physically justified assumptions, ordered by strength, that will simplify the coming calculations.
-
(i) Assumption 1 (weakest): taking the time average of (2.10a ) shows that
$\boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{u}'\rangle = O(M^3)$
, meaning
$\langle \boldsymbol{u}'\rangle$
corresponds to a steady, viscosity-free incompressible mean flow to the order considered here. Under these circumstances, we would expect this mean flow to be vorticity free, implying
$\boldsymbol{\nabla \langle u'\rangle } - (\boldsymbol{\nabla \langle u'\rangle })^T = O(M^2)$
, and consequently, in light of (2.12), that (2.10b
) may be equivalently written as(2.13)This assumption may also be justified by noting that we would expect the flow in a musical instrument to be vorticity free.
\begin{equation} \frac {\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{\nabla \!} p' = \frac {1}{2}\boldsymbol{\nabla }\big (p'^2 - |\boldsymbol{u}'|^2\big ). \end{equation}
-
(ii) Assumption 2a (stronger): on top of assumption 1, by taking a time average of (2.13) we see that the quantity
$\langle p'\rangle - ({1}/{2})\langle p'^2 - |\boldsymbol{u}'|^2\rangle$
is spatially constant, where the constant corresponds to an
$O(M)$
change in the far-field ambient pressure
$p_0$
. By suitably choosing
$p_0$
we may therefore set this change in ambient pressure to be zero, and consequently
$\langle p'\rangle = O(M^2)$
. -
(iii) Assumption 2b (also stronger): on top of assumption 1, we may assume the mean flow
$\langle \boldsymbol{u'}\rangle$
to be identically zero. This is justified by imposing an inlet condition with no mean flow, and noting that no mean flow is induced by the non-zero frequency modes. Thus,
$\langle \boldsymbol{u}'\rangle = \boldsymbol{0}$
. This agrees with our physical intuition that, when musical instruments are played, no flow of significant Mach number is expected. This assumption would therefore be valid for musical instruments, but perhaps not for, for example, thermoacoustic engines, where the acoustics driving a mean flow is an important component of the engine.
Henceforth, we will take all of these assumptions to be true. Defining the coefficient of nonlinearity
$ \beta_0 := ( \unicode{x03B3} + 1)/2$
, and introducing the
$O(M^2)$
quantity
$Q' := (p'^2 - |\boldsymbol{u}'|^2)/2$
(which could arguably be thought of as a Lagrangian of the perturbation), we have
\begin{align} \frac {\partial p'}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}' = \frac {\partial }{\partial t}\big ( \beta_0 {p}^{\prime 2} - Q'\big ), \quad \frac {\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{\nabla \!} p' = \boldsymbol{\nabla } Q'. \end{align}
Note that
$ \beta_0$
has been written with a zero subscript in order to avoid confusion with the spatial-mode enumerator
$\beta$
used later on.
2.3. Fourier modal decomposition
At this point, because we are interested in tonal acoustics, we assume a base frequency
$\omega$
and decompose pressure and velocity into Fourier modes by writing
\begin{align} p'(\boldsymbol{x},t) = \sum _{a = -\infty }^{\infty }P^a(\boldsymbol{x}){\textrm{e}}^{-{\textrm{i}}\kern0.5pt a\omega t}, \quad \boldsymbol{u}'(\boldsymbol{x},t) = \sum _{a = -\infty }^{\infty }\boldsymbol{U}^a(\boldsymbol{x}){\textrm{e}}^{-{\textrm{i}}\kern0.5pt a\omega t}. \end{align}
Here,
$\omega$
is a dimensionless frequency, representing the Helmholtz number. To ensure that
$p'$
and
$\boldsymbol{u}'$
are real, we must have
$P^{-a} = P^{a*}$
and
$\boldsymbol{U}^{-a} = \boldsymbol{U}^{a*}$
, where an asterisk denotes the complex conjugate. Our equations for
$p'$
and
$\boldsymbol{u}'$
(with
$a=0$
modes discounted due to assumptions 2a and 2b above) become
\begin{equation} -{\textrm{i}}\kern0.5pt a\omega P^a + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}^a = -{\textrm{i}}\kern0.5pt a\omega \bigg (- Q^a + \sum _{b \neq 0,a} \beta_0 P^{a - b}P^b \bigg ), \end{equation}
\begin{equation} -{\textrm{i}}\kern0.5pt a\omega \boldsymbol{U}^a + \boldsymbol{\nabla }\! P^a = \boldsymbol{\nabla }\!Q^a, \end{equation}
where
$Q'$
has a Fourier series coefficient given by
\begin{equation} Q^a = \frac {1}{2}\sum _{b \neq 0,a}P^{a - b}P^b - \boldsymbol{U^{a - b}\boldsymbol{\cdot }U^b}. \end{equation}
For brevity, we will now assume
$P^0 \equiv U^0 \equiv 0$
and simply sum over all
$b$
.
2.4. Duct coordinate system
2.4.1. The coordinate system in two dimensions
The duct geometry in two dimensions (left) and three dimensions (right). A Frenet–Serret frame is employed in each case, with centreline
$\boldsymbol{q}(s)$
, tangent
$\boldsymbol{t}(s)$
and normal
$\boldsymbol{n}(s)$
; additionally, in three dimensions we have a binormal
$\boldsymbol{b}(s)$
. The two-dimensional duct has independently varying wall widths
$X_+(s)$
and
$X_-(s)$
, whereas the three-dimensional duct has a single axisymmetrically varying wall-width
$R(s)$
.
As shown in figure 1 (left), the two-dimensional duct is defined by a centreline
$\boldsymbol{q}(s)$
(where
$s$
is the arclength parameterising distance along the duct), a total width at each point
$X(s)$
, and independently varying wall widths
$X_-(s) \lt 0$
and
$X_+(s) \gt 0$
, defined such that
$X(s) = X_+(s) - X_-(s)$
. (The choice to include independently varying wall widths was made early in the formulation of the model, in order to increase the degrees of freedom; in this paper, however, we take the centreline to remain in the centre throughout. Asymmetric wall-width variation would be an interesting area of future study, for example for ducts with cavities on only one side.) The direction along the duct is the longitudinal direction, and the direction perpendicular to this is the transverse direction. The longitudinal direction has coordinate
$s$
, and for the transverse direction we introduce a coordinate
$x \in [X_-(s), X_+(s)]$
. This coordinate basis arises from the Frenet–Serret frame defined by the centreline
$\boldsymbol{q}(s)$
, i.e. the tangent and normal unit vectors
$\boldsymbol{t}$
and
$\boldsymbol{n}$
satisfying
\begin{align} \frac {{\rm d} \boldsymbol{q}}{{\rm d} s} = \boldsymbol{t}, \quad \frac {{\rm d} \boldsymbol{t}}{{\rm d} s} = \kappa \boldsymbol{n}, \quad \frac {{\rm d} \boldsymbol{n}}{{\rm d} s} = -\kappa \boldsymbol{t}, \end{align}
with the scalar
$\kappa (s)$
representing the curvature of the duct. Any point in the duct then has coordinates
$(s,x)$
and position vector
\begin{equation} \boldsymbol{x} = \boldsymbol{q}(s) + x\boldsymbol{n}(s). \end{equation}
A differential
${\rm d} \boldsymbol{x}$
and a metric
${\rm d} \boldsymbol{x}\boldsymbol{\cdot }{\rm d} \boldsymbol{x}$
may then be calculated, from which we may read off the Lamé coefficients and
$(s,x)$
basis vectors
\begin{align} h_s = 1 - \kappa x, \quad h_x = 1, \quad \boldsymbol{e}_s = \boldsymbol{t}, \quad \boldsymbol{e}_x = \boldsymbol{n}. \end{align}
This is all of the machinery necessary to project the governing (2.16). Setting
$\boldsymbol{U}^a = U^a\boldsymbol{e}_s + V^a\boldsymbol{e}_x$
, we get
\begin{gather} -{\textrm{i}}\kern0.5pt a\omega P^a + \frac {1}{h_s}\frac {\partial U^a}{\partial s} + \frac {1}{h_s}\frac {\partial (h_sV^a)}{\partial x} = -{\textrm{i}}\kern0.5pt a\omega \bigg (\sum _b {\unicode{x03B2}}_0 P^{a - b}P^b - Q^a\bigg ), \end{gather}
\begin{gather} -{\textrm{i}}\kern0.5pt a\omega U^a + \frac {1}{h_s}\frac {\partial P^a}{\partial s} = \frac {1}{h_s}\frac {\partial Q^a}{\partial s}, \qquad \qquad \qquad -{\textrm{i}}\kern0.5pt a\omega V^a + \frac {\partial P^a}{\partial x} = \frac {\partial Q^a}{\partial x}. \end{gather}
We now proceed to eliminate the transverse velocities. Since we are neglecting
$O(M^3)$
or smaller terms, we can form an expression for
$V^a$
in terms of the other acoustic variables
\begin{align} V^a = \frac {1}{{\textrm{i}}\kern0.5pt a\omega }\frac {\partial }{\partial x}\left (P^a - Q^a\right ), \quad {\rm with} \quad Q^a = \frac {1}{2}\sum _b\left (\!\!P^{a - b}P^b - U^{a - b}U^b + \frac {\partial _xP^{a - b}\partial _xP^b}{(a - b)b\omega ^2}\right )\!, \end{align}
where
$\partial _x$
denotes
$\partial /\partial x$
and the error in this expression for
$Q^a$
is
$O(M^3)$
. Since our eventual goal is to obtain ODEs in
$s$
for spatial-mode coefficients of
$P^a$
and
$U^a$
, these will be simplest to implement if we have as few
$s$
-derivatives present as possible. Therefore, we also wish to eliminate
$s$
-derivatives from the right-hand sides of our equations, using
$O(M)$
substitutions of the form
\begin{equation} \frac {\partial P^a}{\partial s} = {\textrm{i}}\kern0.5pt a\omega h_sU^a + O(M^2). \end{equation}
Once we have eliminated
$V^a$
and all
$s$
-derivatives at second order, we are left with
\begin{align} \frac {\partial U^a}{\partial s} -& {\textrm{i}}\kern0.5pt a\omega \left [h_s\left (1 + \frac {1}{a^2\omega ^2}\frac {\partial ^2}{\partial x^2}\right ) - \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial x}\right ]P^a \notag\\=& {\textrm{i}}\kern0.5pt a\omega \Bigg \{- {\unicode{x03B2}}_0 h_s\sum _bP^{a - b}P^b + \bigg [h_s\left (1 - \frac {1}{a^2\omega ^2}\frac {\partial ^2}{\partial x^2}\right ) + \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial x}\bigg ]Q^a\Bigg \}, \\[-12pt]\nonumber \end{align}
\begin{align} \frac {\partial P^a}{\partial s} -& {\textrm{i}}\kern0.5pt a\omega h_sU^a \notag\\ &\quad = {\textrm{i}}\kern0.5pt\omega h_s\sum _b\Bigg \{U^{a - b}\left [(a - b) - b\left (1 + \frac {1}{b^2\omega ^2}\frac {\partial ^2}{\partial x^2}\right )\right ]P^b + \frac {\partial _xU^{a - b}\partial _xP^b}{b\omega ^2}\Bigg \}. \end{align}
The other piece of information that needs expression in the coordinate system is the hard-walled boundary condition
$\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nu }_\pm = 0$
at
$x = X_\pm$
, where
$\boldsymbol{\nu }_\pm$
is the normal to each duct wall (distinct from the Frenet–Serret normal
$\boldsymbol{n}$
). The duct walls are defined by the equations
$x - X_\pm (s) = 0$
, so the normals are
$\boldsymbol{\nu }_\pm = \boldsymbol{\nabla }(x - X_\pm ) = \boldsymbol{e}_x - (\partial _sX_\pm /h_s)\boldsymbol{e}_s$
. The boundary condition is then given by the equation
\begin{equation} U^a\frac {\partial X_\pm }{\partial s} = h_sV^a = \frac {h_s}{{\textrm{i}}\kern0.5pt a\omega }\frac {\partial }{\partial x}\left (P^a - Q^a\right ) \quad \text{at}\quad x = X_\pm (s). \end{equation}
2.4.2. The coordinate system in three dimensions
Later on we will be projecting the equations onto a basis of functions representing the modes of a straight duct. The construction of such a basis is much easier and cleaner with separable boundary conditions in three dimensions, so here we drop the asymmetry about the centreline to focus on three-dimensional ducts with circular cross-sections; a constraint most brass instruments adhere to. However, as well as bending in the plane, the duct may now twist out of it, complicating the Frenet–Serret frame.
As shown in figure 1 (right), we begin once more with a centreline
$\boldsymbol{q}(s)$
, but this time have a single radius function
$R(s)$
. The longitudinal coordinate is
$s$
as before, and two more variables span the plane perpendicular to the centreline tangent. The Frenet–Serret frame now has an extra vector and a corresponding extra scalar
\begin{align} \frac {{\rm d} \boldsymbol{q}}{{\rm d} s} = \boldsymbol{t}, \quad \frac {{\rm d} \boldsymbol{t}}{{\rm d} s} = \kappa \boldsymbol{n}, \quad \frac {{\rm d} \boldsymbol{n}}{{\rm d} s} = -\kappa \boldsymbol{t} + \tau \boldsymbol{b}, \quad \frac {{\rm d} \boldsymbol{b}}{{\rm d} s} = -\tau \boldsymbol{n}. \end{align}
Here,
$\boldsymbol{b}$
is the binormal (satisfying
$\boldsymbol{b} = \boldsymbol{t} \times \boldsymbol{n}$
) and
$\tau$
is the torsion. We now introduce polar coordinates
$r \in [0,R(s)]$
and
$\theta \in [0,2 \pi)$
in the transverse plane, but rather than rotating this polar frame in line with the rotation of the basis vectors
$\boldsymbol{n}$
and
$\boldsymbol{b}$
, we leave an extra degree of freedom here by also introducing a phase shift, or twist,
$\theta _0(s)$
, that may vary along the duct (following the work of Germano Reference Germano1982), so that any point
$\boldsymbol{x}$
is given in these coordinates by
\begin{equation} \boldsymbol{x} = \boldsymbol{q}(s) + r\cos (\theta - \theta _0)\boldsymbol{n} + r\sin (\theta - \theta _0)\boldsymbol{b}. \end{equation}
The corresponding metric
${\rm d} \boldsymbol{x}\boldsymbol{\cdot }{\rm d} \boldsymbol{x}$
is then
\begin{align} {\rm d} \boldsymbol{x}\boldsymbol{\cdot }{\rm d} \boldsymbol{x} = {\rm d} s^2\big [1 - \kappa r\cos (\theta - \theta _0) + r^2(\theta _0' - \tau )^2\big ] + {\rm d} r^2 + r^2\,{\rm d} \theta ^2 - 2\,{\rm d} s\,{\rm d} \theta r^2(\theta _0' - \tau ). \end{align}
An orthogonal coordinate system requires that there be no cross-term differentials in the metric: to achieve this, following Germano (Reference Germano1982), we utilise the extra degree of freedom and take
$\theta _0' = \tau$
. The Lamé coefficients and
$(s,r,\theta )$
basis vectors can then be written as
\begin{align} h_s &= 1 - \kappa r\cos (\theta -\theta _0), \quad h_r = 1, \quad h_\theta = r, \quad \boldsymbol{e}_s = \boldsymbol{t}, \end{align}
\begin{align} \boldsymbol{e}_r &= \boldsymbol{n}\cos (\theta -\theta _0) + \boldsymbol{b}\sin (\theta -\theta _0), \quad \boldsymbol{e}_\theta = -\boldsymbol{n}\sin (\theta -\theta _0) + \boldsymbol{b}\cos (\theta -\theta _0). \end{align}
We decompose the velocity in the same way as before, but now with an extra coordinate, i.e.
$\boldsymbol{U}^a = U^a\boldsymbol{e}_s + V^a\boldsymbol{e}_r + W^a\boldsymbol{e}_\theta$
. The mass and momentum (2.16) become
\begin{gather} -{\textrm{i}}\kern0.5pt a\omega P^a + \frac {1}{h_s}\frac {\partial U^a}{\partial s} + \frac {1}{rh_s}\frac {\partial (rh_sV^a)}{\partial r} + \frac {1}{rh_s}\frac {\partial (h_sW^a)}{\partial \theta } = {\textrm{i}}\kern0.5pt a\omega \bigg (\!- \unicode{x03B2}_0\sum _bP^{a - b}P^b + Q^a\bigg ), \end{gather}
\begin{gather} -{\textrm{i}}\kern0.5pt a\omega U^a + \frac {1}{h_s}\frac {\partial P^a}{\partial s} = \frac {1}{h_s}\frac {\partial Q^a}{\partial s}, \quad -{\textrm{i}}\kern0.5pt a\omega V^a + \frac {\partial P^a}{\partial r} = \frac {\partial Q^a}{\partial r}, \quad -{\textrm{i}}\kern0.5pt a\omega W^a + \frac {1}{r}\frac {\partial P^a}{\partial \theta } = \frac {1}{r}\frac {\partial Q^a}{\partial \theta }. \end{gather}
Proceeding as before, we find expressions for
$V^a$
and
$W^a$
in terms of
$P^a$
and
$U^a$
\begin{align} V^a = \frac {1}{{\textrm{i}}\kern0.5pt a\omega }\frac {\partial }{\partial r}(P^a - Q^a), \quad W^a = \frac {1}{{\textrm{i}}\kern0.5pt a\omega r}\frac {\partial }{\partial \theta }(P^a - Q^a), \end{align}
with
\begin{equation} Q^a = \frac {1}{2}\sum _b\left (P^{a - b}P^b - U^{a - b}U^b + \frac {\boldsymbol{\nabla }_{{t}}P^{a - b}\boldsymbol{\cdot }\boldsymbol{\nabla }_{{t}}P^b}{(a - b)b\omega ^2}\right )\!, \end{equation}
where we have introduced the transverse gradient
$\boldsymbol{\nabla }_{{t}} = \boldsymbol{e_r}\partial _r + \tfrac {1}{r}\boldsymbol{e_\theta }\partial _\theta$
which acts only in the
$r,\theta$
plane normal to the centreline. We also define, for any function
$f$
, the Frenet–Serret normal derivative (FSND)
\begin{equation} \frac {\partial f}{\partial \boldsymbol{n}} = \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\! f = \cos (\theta -\theta _0)\frac {\partial f}{\partial r} - \frac {\sin (\theta -\theta _0)}{r}\frac {\partial f}{\partial \theta }. \end{equation}
This is a useful quantity, since it allows us to commute the transverse gradient with the longitudinal scale factor
$h_s$
, as in
\begin{align} &\boldsymbol{\nabla }_{{t}}\boldsymbol{\cdot }\left (h_s\boldsymbol{\nabla }_{{t}}f\right ) = \big (h_s\boldsymbol{\nabla }_{{t}}^2 - \kappa \partial _{\boldsymbol{n}}\big )f \quad {\rm and} \quad {\nabla }_{{t}}^2(h_sf) = \big (h_s{\nabla }_{{t}}^2 - 2\kappa \partial _{\boldsymbol{n}}\big )f. \end{align}
This contrasts with the two-dimensional derivation, where both the transverse gradient and the FSND collapse onto the
$x$
-derivative (the former as a vector in the
$\boldsymbol{e}_x$
direction, the latter as a scalar). Once we have eliminated the
$s$
-derivatives as before, we then have
\begin{align} &\frac {\partial U^a}{\partial s} - {\textrm{i}}\kern0.5pt a\omega \left [h_s\! \left (1 + \frac {{\nabla }_{{t}}^2}{a^2\omega ^2}\right ) - \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial \boldsymbol{n}}\right ]P^a \nonumber \\ &\qquad = {\textrm{i}}\kern0.5pt a\omega \Bigg \{\!- \beta_0 h_s\sum _bP^{a - b}P^b + \bigg [h_s\! \left (1 - \frac {{\nabla }_{{t}}^2}{a^2\omega ^2}\right ) + \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial \boldsymbol{n}}\bigg ]Q^a\Bigg \}, \end{align}
\begin{align} &\frac {\partial P^a}{\partial s} - {\textrm{i}}\kern0.5pt a\omega h_sU^a = {\textrm{i}}\kern0.5pt\omega h_s\sum _b\Bigg \{U^{a - b}\!\left [(a - b) - b\!\left (1 + \frac {{\nabla }_{{t}}^2}{b^2\omega ^2}\right )\right ]P^b + \frac {\boldsymbol{\nabla }_{{t}}U^{a - b}\boldsymbol{\cdot }\boldsymbol{\nabla }_{{t}}P^b}{b\omega ^2}\Bigg \}. \end{align}
Finally, the hard-walled boundary condition in three dimensions requires the normal
$\boldsymbol{\nu }$
, which we calculate similarly to earlier as
$\boldsymbol{\nu } = \boldsymbol{\nabla }(r - R)|_{r = R} = (\boldsymbol{e}_r - R'\boldsymbol{e}_s/h_s)|_{r = R}$
. Plugging this into
$\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nu }|_{r = R} = 0$
gives
\begin{equation} R'U^a|_{r = R} = (h_sV^a)|_{r = R} = \left [\frac {h_s}{{\textrm{i}}\kern0.5pt a\omega }\frac {\partial }{\partial r}\left (P^a - Q^a\right )\right ]\bigg |_{r = R}. \end{equation}
2.5. Spatial modal decomposition
At this point, we expand each of
$P^a$
and
$U^a$
in terms of a basis of straight-duct modes. We will ultimately be solving for the coefficients of the series expansions, which will only depend on
$s$
with all transverse variation being contained within the basis functions.
2.5.1. The spatial modes in two dimensions
The variables
$P^a$
and
$U^a$
are expanded as
\begin{align} &P^a = \sum _{\alpha = 0}^\infty P_\alpha ^a(s)\psi _\alpha (s,x), \quad U^a = \sum _{\alpha = 0}^\infty U_\alpha ^a(s)\psi _\alpha (s,x), \end{align}
with every
$\psi _\alpha$
satisfying
-
(i) Helmholtz’s equation with eigenvalue
$\lambda _\alpha$
in two dimensions (scaling out the eigenvalue’s
$X$
-dependence for convenience)(2.38)
\begin{equation} \frac {\partial ^2\psi _\alpha }{\partial x^2} + \frac {\lambda _\alpha ^2}{X^2}\psi _\alpha = 0; \end{equation}
-
(ii) a normalisation condition
(2.39)
\begin{equation} \langle \psi _\alpha ,\psi _\beta \rangle = \delta _{\alpha \beta }, \quad {\rm where } \quad \langle \psi _\alpha ,\psi _\beta \rangle := \int _{X_-}^{X_+}\psi _\alpha \psi _\beta \,{\rm d} x; \end{equation}
-
(iii) and a Neumann condition on the duct walls (which ensures that the requisite Sturm–Liouville properties hold, while maintaining consistency with the no-penetration condition)
(2.40)
\begin{equation} \frac {\partial \psi _\alpha }{\partial x}\bigg |_{x = X_\pm } = 0. \end{equation}
Solutions to the Helmholtz equation are quantised by the boundary conditions, giving solutions
\begin{align} \psi _\alpha = \frac {C_\alpha }{\sqrt {X}} \cos\! \left [\frac {\lambda _\alpha (x - X_-)}{X}\right ], \quad \lambda _\alpha = \alpha \pi, \end{align}
for
$\alpha \in \mathbb{N}_0$
, with
$C_\alpha$
to be
$\sqrt {2 - \delta _{\alpha 0}}$
ensuring orthonormality.
2.5.2. The spatial modes in three dimensions
In three dimensions, the expansions are similarly written as
\begin{align} P^a = \sum _{\alpha = 0}^\infty P_\alpha ^a(s)\psi _\alpha (s,r,\theta ), \quad U^a = \sum _{\alpha = 0}^\infty U_\alpha ^a(s)\psi _\alpha (s,r,\theta ), \end{align}
with every
$\psi _\alpha$
satisfying the Helmholtz equation in the circular cross-section with Neumann boundary conditions. This gives
\begin{equation} \psi _\alpha = \frac {C_\alpha }{\sqrt { \pi}R} {{J}}_p(\lambda _\alpha r/R)\cos\! \left (p_\alpha (\theta -\theta _0) - \frac {\xi _\alpha \pi}{2}\right )\!, \end{equation}
where
$C_\alpha$
ensures orthonormalisation,
$\lambda _\alpha$
is the eigenvalue,
$p_\alpha \in \mathbb{N}_0$
gives the azimuthal order and
$\xi _\alpha \in \{0,1\}$
effectively provides both sine and cosine solutions. Defining
$j_{pq}'$
to be the
$q^{{}}$
th zero of
${{J}}_p'(x)$
, we have
$\lambda _\alpha = j_{p_\alpha q_\alpha }'$
for
$q_\alpha \in \mathbb{N}_0$
independent of
$\xi _\alpha$
. Hence, modes are given by a bijection between
$\alpha \in \mathbb{N}_0$
and
$(p_\alpha ,q_\alpha ,\xi _\alpha )\in \mathbb{N}_0\times \mathbb{N}_0\times \{0,1\}$
, although modes that are identically zero are skipped (such as
$(0,0,1)$
,
$(2,0,0)$
,
$(2,0,1)$
,
$(4,0,0)$
, etc.). For orthonormality, we have
\begin{align} \langle \psi _\alpha ,\psi _\beta \rangle = \!\!\int _{0}^{R}\!\!\!\!\int _0^{2 \pi}\!\!\!\!\psi _\alpha \psi _\beta \,r{\rm d} r{\rm d}\theta \quad \Rightarrow \quad C_\alpha = \begin{cases} \big |{{J}}_{0}(j_{0 q_\alpha }')\big |^{-1},&p_\alpha = 0,\\[5pt] \left (\sqrt {\dfrac {1}{2}\left [1 - \frac {{p_\alpha }^2}{{j_{p_\alpha q_\alpha }'}^2}\right ]} \,\big |{{J}}_{p_\alpha }(j_{p_\alpha q_\alpha }')\big |\right )^{-1}\!,&p_\alpha \neq 0. \end{cases} \end{align}
2.6. Spatial projection and notation
Now that we have defined the spatial basis functions, we follow McTavish & Brambley (Reference McTavish and Brambley2019) in projecting the governing (2.24, 2.25) or (2.35, 2.36) onto the spatial modal basis (2.41) or (2.43) to form an infinite set of ODEs for the coefficients
$U^a_\alpha$
and
$P^a_\alpha$
. This is achieved by multiplying each of the governing equations by the spatial basis functions
$\psi _\beta$
and integrating across a duct cross-section, using the orthogonality of the spatial basis functions in the process. The result is simplified greatly using a compact notation, which we now introduce.
For each
$a\in \mathbb{N}$
, let
$\boldsymbol{p}^a$
denote the vector of coefficients
$P^a_\alpha$
for
$\alpha \in \mathbb{N}_0$
, and similarly let
$\boldsymbol{u}^a$
denote the vector of coefficients
$U^a_\alpha$
. We use Roman letters as superscripts for the temporal Fourier modal decomposition, and Greek letters as subscripts for the spatial modal decomposition. For each of these vectors, matrices
$ \unicode{x1D648}$
with coefficients
$ \unicode{x1D648}_{\alpha \beta }$
act in the normal way
\begin{equation} \big ( \unicode{x1D648}\boldsymbol{y}\big )_\alpha = \sum _{\beta =0}^\infty \unicode{x1D648}_{\alpha \beta }y_\beta . \end{equation}
For the weakly nonlinear terms, we define a quadratic operator
$\mathcal{M}$
with coefficients
$\mathcal{M}_{\alpha \beta \gamma }$
, which acts on vectors
$\boldsymbol{y}$
and
$\boldsymbol{z}$
as
\begin{equation} \big (\mathcal{M}\langle \boldsymbol{y},\boldsymbol{z}\rangle \big )_\alpha = \sum _{\beta =0}^\infty \sum _{\gamma =0}^\infty \mathcal{M}_{\alpha \beta \gamma }y_\beta z_\gamma . \end{equation}
We also use the operator shorthand
\begin{equation} \big (\mathcal{M} + \mathcal{N}\langle \unicode{x1D644}, \unicode{x1D63C}\rangle \big )\langle \boldsymbol{y},\boldsymbol{z}\rangle = \mathcal{M}\langle \boldsymbol{y},\boldsymbol{z}\rangle + \mathcal{N}\langle \boldsymbol{y}, \unicode{x1D63C}\boldsymbol{z}\rangle . \end{equation}
Unlike in McTavish & Brambley (Reference McTavish and Brambley2019), here, we choose a decomposition into matrix and quadratic operators that depends only on the modal basis and not on the variable
$s$
or any physical parameters. This allows for a unified approach to both the two- and three-dimensional cases listed above, and, since physical quantities such as curvature
$\kappa$
occur explicitly, the effects of the physical quantities can be isolated and better understood. The algebra behind the manipulations needed to get to the governing equations is given in Appendix A. In two dimensions, this results in
\begin{gather} \begin{aligned} \Bigg [\frac {{\rm d} }{{\rm d} s} &+ \frac {X'(s)}{2X(s)} \unicode{x1D652} + \frac {X_-'(s)}{X(s)}\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}\Bigg ]\boldsymbol{u}^a(s) \\ &- {\textrm{i}}\kern0.5pt a\omega \left [\left ( \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X(s)^2}\right )\big ((1 - \kappa (s) X_-(s)) \unicode{x1D644} - \kappa (s) X(s) \unicode{x1D63C}\big ) - \frac {\kappa (s)\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}}{a^2\omega ^2X(s)}\right ]\boldsymbol{p}^a(s) \\ =& \frac {{\textrm{i}}\kern0.5pt a\omega }{\sqrt {X(s)}}\sum _b\left \{\vphantom{\frac {\substack{ \Bigg ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X(s)^2}\Bigg )\big ((1 - \kappa (s) X_-(s))\mathcal{I}^\lambda - \kappa (s) X(s)\mathcal{A}^\lambda \big ) + \dfrac {\kappa (s)\kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda }{a^2\omega ^2X(s)}}}{2(a - b)b\omega ^2X(s)^2}}- \beta_0\big ((1 - \kappa (s) X_-(s))\mathcal{I} - \kappa X(s)\mathcal{A}\big )\langle \boldsymbol{p}^{a - b}(s),\boldsymbol{p}^b(s)\rangle\right. \\& +\Bigg [\Bigg ( \unicode{x1D644} + \frac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X(s)^2}\Bigg )\big ((1 - \kappa (s) X_-(s))\mathcal{I} - \kappa (s) X(s)\mathcal{A}\big )\\& \quad + \frac {\kappa (s)\kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}}{a^2\omega ^2X(s)}\Bigg ]\frac {\langle \boldsymbol{p}^{a - b}(s),\boldsymbol{p}^b(s)\rangle - \langle \boldsymbol{u}^{a - b}(s),\boldsymbol{u}^b(s)\rangle }{2} \\& + \left [\frac {{ \Bigg ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X(s)^2}\Bigg )\big ((1 - \kappa (s) X_-(s))\mathcal{I}^\lambda - \kappa (s) X(s)\mathcal{A}^\lambda \big ) + \dfrac {\kappa (s)\kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda }{a^2\omega ^2X(s)}}}{2(a - b)b\omega ^2X(s)^2} \right ]\\ &\left.\quad\times\langle \boldsymbol{p}^{a - b}(s),\boldsymbol{p}^b(s)\rangle \vphantom{\frac {\substack{ \Bigg ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X(s)^2}\Bigg )\big ((1 - \kappa (s) X_-(s))\mathcal{I}^\lambda - \kappa (s) X(s)\mathcal{A}^\lambda \big ) + \dfrac {\kappa (s)\kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda }{a^2\omega ^2X(s)}}}{2(a - b)b\omega ^2X(s)^2}}\right \}, \end{aligned} \end{gather}
\begin{gather} \begin{aligned} \Bigg [\frac {{\rm d} }{{\rm d} s} - &\frac {X'(s)}{2X(s)} \unicode{x1D652}^{{\,T}} - \frac {X_-'(s)}{X(s)}\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}^{{T}}\Bigg ]\boldsymbol{p}^a(s) - {\textrm{i}}\kern0.5pt a\omega \big [(1 - \kappa (s) X_-(s)) \unicode{x1D644} - \kappa (s) X(s) \unicode{x1D63C}\big ]\boldsymbol{u}^a(s) \\[-9pt]\\ =& \frac {1}{\sqrt {X(s)}}\sum _b\Bigg \{\bigg (\frac {X_+'(s)}{X(s)}\overline {\mathcal{W}}^+ - \frac {X_-'(s)}{X(s)}\overline {\mathcal{W}}^-\bigg )\langle \boldsymbol{u}^{a - b}(s),\boldsymbol{u}^b(s)\rangle \\[-9pt]\\& + {\textrm{i}}\kern0.5pt\omega \Bigg [\big ((1 - \kappa (s) X_-(s))\mathcal{I} - \kappa (s) X(s)\mathcal{A}\big ) \left \langle \unicode{x1D644},(a - b) \unicode{x1D644} - b\left ( \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{b^2\omega ^2X(s)^2}\right )\right \rangle \\[-9pt]\\&\quad + \frac {(1 - \kappa (s) X_-(s))\mathcal{I}^\lambda - \kappa (s) X(s)\mathcal{A}^\lambda }{b\omega ^2X(s)^2}\Bigg ]\langle \boldsymbol{u}^{a - b}(s),\boldsymbol{p}^b(s)\rangle \Bigg \}, \end{aligned} \end{gather}
where the matrices and tensors are defined in Appendix A.1 and are independent of
$s$
.
Likewise, in three dimensions, using the derivation in Appendix A.2, we arrive at
\begin{align} &\left (\frac {{\rm d} }{{\rm d} s} + \frac {R'(s)}{R(s)} \unicode{x1D652} + \tau (s) \unicode{x1D643}\right )\boldsymbol{u}^a(s) \notag\\ \notag &\qquad \quad - {\textrm{i}}\kern0.5pt a\omega \left [\left ( \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R(s)^2}\right )\left ( \unicode{x1D644} - \kappa (s) R(s) \unicode{x1D63C}\right ) - \frac {\kappa (s)\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}}{a^2\omega ^2R(s)}\right ]\boldsymbol{p}^a(s) \\ \notag &\quad = \frac {{\textrm{i}}\kern0.5pt a\omega }{\sqrt { \pi}R(s)}\sum _b\Bigg \{ \Bigg [\Bigg ( \unicode{x1D644} + \frac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R(s)^2}\Bigg )\left (\mathcal{I} - \kappa (s) R(s)\mathcal{A}\right ) + \frac {\kappa (s)\kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}}{a^2\omega ^2 R(s)}\Bigg ] \\\notag & \quad\quad\quad\quad\quad \times \frac {\langle \boldsymbol{p}^{a - b}(s),\boldsymbol{p}^b(s)\rangle - \langle \boldsymbol{u}^{a - b}(s),\boldsymbol{u}^b(s)\rangle }{2} \\\notag &\quad\quad\quad\quad + \Bigg [\frac {\bigg ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R(s)^2}\bigg )\left (\mathcal{I}^\lambda - \kappa (s) R(s)\mathcal{A}^\lambda \right ) + \dfrac {\kappa (s)\kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda }{a^2\omega ^2R(s)}}{2(a - b)b\omega ^2R(s)^2} \\&\quad\quad\quad\quad - \beta_0\left (\mathcal{I} - \kappa (s) R(s)\mathcal{A}\right )\Bigg ] \langle \boldsymbol{p}^{a - b}(s),\boldsymbol{p}^b(s)\rangle \Bigg \}, \end{align}
\begin{align} &\left (\frac {{\rm d} }{{\rm d} s} - \frac {R'(s)}{R(s)} \unicode{x1D652}^{{\,T}} - \tau (s) \unicode{x1D643}^{{T}}\right )\boldsymbol{p}^a(s) - {\textrm{i}}\kern0.5pt a\omega \big ( \unicode{x1D644} - \kappa (s) R(s) \unicode{x1D63C}\big )\boldsymbol{u}^a(s) \notag\\ \notag &\qquad = \frac {1}{\sqrt { \pi}R(s)}\sum _b\Bigg \{\frac {R'(s)}{R(s)}\overline {\mathcal{W}}\langle \boldsymbol{u}^{a - b}(s),\boldsymbol{u}^b(s)\rangle + {\textrm{i}}\kern0.5pt\omega \Bigg [\left (\mathcal{I} - \kappa (s) R(s)\mathcal{A}\right ) \\&\qquad \quad \times \bigg \langle \unicode{x1D644},(a - b) \unicode{x1D644} - b\bigg ( \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{b^2\omega ^2R(s)^2}\bigg )\bigg \rangle + \frac {\mathcal{I}^\lambda - \kappa (s) R(s)\mathcal{A}^\lambda }{b\omega ^2R(s)^2}\Bigg ]\langle \boldsymbol{u}^{a - b}(s),\boldsymbol{p}^b(s)\rangle \Bigg \}, \end{align}
where again the matrices and tensors are defined in Appendix A.2 and are independent of
$s$
.
The equations have a structure common to both two and three dimensions, and can be written compactly (as detailed in Appendix B) as
\begin{equation} \frac {{\rm d} }{{\rm d} s}\!\begin{pmatrix} \boldsymbol{u}^a \\ \boldsymbol{p}^a \end{pmatrix} = \unicode{x1D647}^a\! \begin{pmatrix} \boldsymbol{u}^a \\ \boldsymbol{p}^a \end{pmatrix} + \sum _b\mathcal{N}^{\,ab}\! \left \langle \begin{pmatrix} \boldsymbol{u}^{a - b} \\ \boldsymbol{p}^{a - b} \end{pmatrix},\begin{pmatrix} \boldsymbol{u}^b \\ \boldsymbol{p}^b \end{pmatrix}\right \rangle \!, \end{equation}
in terms of square blocks
$ \unicode{x1D647}^a =: \begin{pmatrix} \unicode{x1D647}_1^a & \unicode{x1D647}_2^a \\ \unicode{x1D647}_3^a & \unicode{x1D647}_4^a \end{pmatrix}$
, and ‘cubic blocks’
\begin{equation} \mathcal{N}^{\,ab}\! \left \langle \begin{pmatrix} \boldsymbol{u}^{a - b} \\ \boldsymbol{p}^{a - b} \end{pmatrix},\begin{pmatrix} \boldsymbol{u}^b \\ \boldsymbol{p}^b \end{pmatrix}\right \rangle = \begin{pmatrix} \mathcal{N}^{\,ab}_1\langle \boldsymbol{u}^{a-b}, \boldsymbol{u}^b\rangle + \mathcal{N}^{\,ab}_6\langle \boldsymbol{p}^{a-b}, \boldsymbol{p}^b\rangle \\[5pt] \mathcal{N}^{\,ab}_3\langle \boldsymbol{u}^{a-b}, \boldsymbol{u}^b\rangle + \mathcal{N}^{\,ab}_7\langle \boldsymbol{u}^{a-b}, \boldsymbol{p}^b\rangle \end{pmatrix}\!, \end{equation}
where the equivalent terms
$\mathcal{N}^{\,ab}_2 \equiv \mathcal{N}^{\,ab}_4 \equiv \mathcal{N}^{\,ab}_5 \equiv \mathcal{N}^{\,ab}_8 \equiv 0$
. We recall that
$P_\alpha ^{-a} = P_\alpha ^{a*}$
and that
$P_\alpha ^0=0$
by our earlier assumptions of real variables and vanishing time averages. Note that, as expected, the weakly nonlinear terms will mix the effects of different frequencies. We will now work in terms of these quantities, so that everything we do will apply in both two and three dimensions.
In summary, the governing equations are (2.48) in two dimensions and (2.49) in three dimensions, and can be combined and written as (2.50). Greater numerical efficiency has been achieved through the transferral of
$s$
-dependence from matrices to the scalars multiplying them (since this allows the matrices to be pre-calculated rather than requiring recalculation at each step). Upper (Roman) indices are Fourier series modenumbers while lower (Greek) indices are duct modenumbers.
3. Solution of the governing equations
Here, we solve the above governing equations using a multi-modal method, inspired by Félix & Pagneux (Reference Félix and Pagneux2001) and McTavish & Brambley (Reference McTavish and Brambley2019). We first truncate the number of spatial modes to
$\alpha _{\textrm{max}}$
and the number of temporal modes to
$a_{\textrm{max}}$
, before introducing the concept of an admittance which relates velocities to pressures.
3.1. Admittance
The admittance,
$Y$
, and its inverse the impedance,
$Z$
, characterise the relationship between acoustic pressure and acoustic velocity. The admittance and impedance are defined here as
\begin{align} &\boldsymbol{u}^a = \unicode{x1D654}^{\,a}\boldsymbol{p}^a + \sum _b\mathcal{Y}^{\,ab}\langle \boldsymbol{p}^{a - b},\boldsymbol{p}^b\rangle , \quad \boldsymbol{p}^a = \unicode{x1D655}^{\,a}\boldsymbol{u}^a + \sum _b\mathcal{Z}^{ab}\langle \boldsymbol{u}^{a - b},\boldsymbol{u}^b\rangle . \end{align}
The first of these relations may be used to eliminate the velocity from the governing (2.50), which results in a Riccati-style equation for the
$s$
-evolution of the admittance
\begin{gather} \frac {{\rm d} \unicode{x1D654}^{\,a}}{{\rm d} s} = - \unicode{x1D654}^{\,a} \unicode{x1D647}_3^a \unicode{x1D654}^{\,a} + \unicode{x1D647}_1^a \unicode{x1D654}^{\,a} - \unicode{x1D654}^{\,a} \unicode{x1D647}_4^a + \unicode{x1D647}_2^a, \end{gather}
\begin{align} \frac {{\rm d} \mathcal{Y}^{\,ab}}{{\rm d} s} = &- \unicode{x1D654}^{\,a}\mathcal{N}_3^{ab}\langle \unicode{x1D654}^{\,a - b}, \unicode{x1D654}^{\,b}\rangle - \unicode{x1D654}^{\,a}\mathcal{N}_7^{ab}\langle \unicode{x1D654}^{\,a - b}, \unicode{x1D644}\,\rangle + \mathcal{N}_1^{ab}\langle \unicode{x1D654}^{\,a - b}, \unicode{x1D654}^{\,b}\rangle + \mathcal{N}_6^{ab} \nonumber\\ &- \mathcal{Y}^{\,ab}\big [- \unicode{x1D647}_1^a + \unicode{x1D654}^{\,a} \unicode{x1D647}_3^a + \big \langle \unicode{x1D647}_4^{a - b} + \unicode{x1D647}_3^{a - b} \unicode{x1D654}^{\,a - b}, \unicode{x1D644}\big \rangle + \big \langle \unicode{x1D644}, \unicode{x1D647}_4^b + \unicode{x1D647}_3^{b} \unicode{x1D654}^{\,b}\big \rangle \big ].\end{align}
These equations will be solved to find the admittance, which encodes the acoustic properties of the duct. Once the admittance is known, the acoustic pressure may be solved for by substituting the admittance definition (3.1) back into the governing (2.50), giving
\begin{align} \begin{aligned} \frac {{\rm d}\boldsymbol{p}^a}{{\rm d} s} = \left ( \unicode{x1D647}_3^a \unicode{x1D654}^{\,a} + \unicode{x1D647}_4^a\right )\boldsymbol{p}^a + \sum _b\big ( \unicode{x1D647}_3^a\mathcal{Y}^{\,ab} + \mathcal{N}_3^{ab}\langle \unicode{x1D654}^{\,a - b}, \unicode{x1D654}^{\,b}\rangle + \mathcal{N}_7^{ab}\langle \unicode{x1D654}^{\,a - b}, \unicode{x1D644}\,\rangle \big )\langle \boldsymbol{p}^{a - b},\boldsymbol{p}^b\rangle , \end{aligned} \end{align}
with the velocities then given by the known admittance and pressure from the admittance definition (3.1). In order to solve the first-order (3.2), a known value of the admittance must be given as a boundary condition at some point in the duct, for example, at the outlet. We next investigate special values of the admittance that might be used as boundary conditions for (3.2).
3.2. Invariant admittances
We consider admittances that solve (3.2) and are constant in
$s$
, therefore having a vanishing
$s$
-derivative. Physically, these invariant solutions represent the admittances of ducts for which no
$s$
position is distinguishable from any other. Such a duct would have a constant radius, curvature and torsion, leaving the possibilities of: (i) an infinite straight duct; (ii) an annulus (two dimensions) or torus (three dimensions); and (iii) an infinite helical duct. Because of the quadratic term in the evolution equation for
$ \unicode{x1D654}^{\,a}$
, there is more than one constant solution for a given geometry. Here, a ‘
$+$
’ solution will be constructed from the set of eigenvalues that represent waves either decaying or propagating in the positive
$s$
-direction, and a ‘
$-$
’ solution will be constructed from the eigenvalues representing growth in the positive
$s$
-direction or propagation in the negative
$s$
-direction. These are respectively the positive and negative characteristic admittances. We denote by
$(\boldsymbol{u}^{a\pm },\boldsymbol{p}^{a\pm })$
disturbances propagating in only the positive
$+$
or negative
$-$
directions, and define
\begin{align} \boldsymbol{u}^{a\pm } &= \unicode{x1D654}^{\,a\pm}\boldsymbol{p}^{a\pm } + \sum _b\mathcal{Y}^{\,ab\pm }\langle \boldsymbol{p}^{(a - b)\pm },\boldsymbol{p}^{b\pm }\rangle , \end{align}
\begin{align} \boldsymbol{p}^{a\pm } &= \unicode{x1D655}^{\,a\pm }\boldsymbol{u}^{a\pm } + \sum _b\mathcal{Z}^{ab\pm }\langle \boldsymbol{u}^{(a - b)\pm },\boldsymbol{u}^{b\pm }\rangle . \end{align}
The linear characteristic impedances
$ \unicode{x1D655}^{\,a\pm }$
are the inverses of
$ \unicode{x1D654}^{\,a\pm}$
; the nonlinear characteristic impedances are found by the weakly nonlinear inversion rule
$\mathcal{Z}^{ab\pm } = - \unicode{x1D655}^{\,a\pm }\mathcal{Y}^{\,ab\pm }\langle \unicode{x1D655}^{(a - b)\pm }, \unicode{x1D655}^{\,b\pm }\rangle$
.
3.2.1. Straight-duct characteristic admittances
This computation will be done in the general case in order to lay the groundwork for the curvature/torsion cases, where the diagonalisation is less trivial.
For a straight duct, the matrix
$ \unicode{x1D647}^a$
is constant, and given by setting
$\kappa \equiv \tau \equiv X' \equiv R' \equiv 0$
\begin{equation} \overline { \unicode{x1D647}}^a = \begin{pmatrix} \boldsymbol{\mathit{0}} &\overline { \unicode{x1D647}}_2^a \\[5pt] \overline { \unicode{x1D647}}_{\,3}^{\,a} &\boldsymbol{\mathit{0}} \end{pmatrix}\!, \end{equation}
where
$\overline { \unicode{x1D647}}_2^a$
is diagonal and
$\overline { \unicode{x1D647}}_{\,3}^{\,a}$
is proportional to the identity. We seek a solution to
\begin{equation} \frac {{\rm d}}{{\rm d} s}\begin{pmatrix} \boldsymbol{u}^a \\ \boldsymbol{p}^a \end{pmatrix} = \overline { \unicode{x1D647}}^a\begin{pmatrix} \boldsymbol{u}^a \\ \boldsymbol{p}^a \end{pmatrix}\!, \end{equation}
which takes the form
$\boldsymbol{p}^a = \overline {\boldsymbol{c}}^a\exp (\overline {\gamma }^a s)$
, with
$\overline {\gamma }^a$
an eigenvalue and
$\overline {\boldsymbol{c}}^a$
the corresponding eigenvector of
$\overline { \unicode{x1D647}}^a$
. We use the fact that for invertible matrices
$ \unicode{x1D63C}$
and
$ \unicode{x1D63F}$
\begin{equation} \det \begin{pmatrix} \unicode{x1D63C} & \unicode{x1D63D} \\ \unicode{x1D63E} & \unicode{x1D63F} \end{pmatrix} = \det \unicode{x1D63C}\det \big ( \unicode{x1D63F} - \unicode{x1D63E} \unicode{x1D63C}^{-1} \unicode{x1D63D}\big ) = \det \unicode{x1D63F}\det \big ( \unicode{x1D63C} - \unicode{x1D63D} \unicode{x1D63F}^{\,-1} \unicode{x1D63E}\big ). \end{equation}
With this in mind, we know that the characteristic equation reads
\begin{equation} 0 = \det\! \left [\overline { \unicode{x1D647}}^a - \overline {\gamma }^a\begin{pmatrix} \unicode{x1D644} & 0 \\ 0 & \unicode{x1D644} \end{pmatrix}\right ] = \det \big (\big (\overline {\gamma }^a\big )^2 \unicode{x1D644} - \overline { \unicode{x1D647}}_{\,3}^{\,a}\overline { \unicode{x1D647}}_2^a\big ). \end{equation}
Because
$\overline { \unicode{x1D647}}_{\,3}^{\,a}\overline { \unicode{x1D647}}_2^a$
is diagonal, we may directly read off the eigenvalues. Since we have truncated to
$\alpha _{\textrm{max}}$
spatial modes, and therefore both
$\overline { \unicode{x1D647}}_2^a$
and
$\overline { \unicode{x1D647}}_{\,3}^{\,a}$
are
$\alpha _{\textrm{max}}\times \alpha _{\textrm{max}}$
, we see that there are precisely
$\alpha _{\textrm{max}}$
solutions for
${\overline {\gamma }^a}^2$
, of which some are positive and some negative (depending on the frequency
$\omega$
, since this is determining whether the modes are cut-on or cutoff). Thus, there will be a set of distinct eigenvalues
$\{\overline {\gamma }_{\,\alpha}^{\,a}\}_{\alpha = 0}^{\alpha _{\textrm{max}}}$
that are exclusively in
$\mathbb{R}^-$
or
${\textrm{i}}\kern0.5pt\mathbb{R}^+$
, corresponding to forward-decaying or forward-propagating modes, respectively. These will be used to build the positive characteristic admittance; each will meanwhile have a mirror image in either
$\mathbb{R}^{+}$
or
${\textrm{i}}\kern0.5pt\mathbb{R}^{-}$
, used to construct the negative characteristic admittance.
We may now partition the eigenvectors into two sets, those corresponding to
$\{\overline {\gamma }_{\,\alpha}^{\,a}\}_{\alpha = 0}^{\alpha _{\textrm{max}}}$
and those corresponding to
$\{-\overline {\gamma }_{\,\alpha}^{\,a}\}_{\alpha = 0}^{\alpha _{\textrm{max}}}$
, denoted
$\overline {\boldsymbol{c}}^{a\pm }$
. Splitting them (as we have split
$\overline { \unicode{x1D647}}^a$
) into upper and lower vectors
$(\overline {\boldsymbol{c}}_{\alpha ,1}^{a\pm },\overline {\boldsymbol{c}}_{\alpha ,2}^{a\pm })$
, and eliminating the upper vectors, we get
\begin{equation} \overline { \unicode{x1D647}}_{\,3}^{\,a}\overline { \unicode{x1D647}}_2^a\overline {\boldsymbol{c}}_{\alpha ,2}^{a\pm } = {\overline {\gamma }_{\,\alpha}^{\,a}}^2\overline {\boldsymbol{c}}_{\alpha ,2}^{a\pm }. \end{equation}
Because
$\overline {\boldsymbol{c}}_{\alpha ,2}^{a+}$
and
$\overline {\boldsymbol{c}}_{\alpha ,2}^{a-}$
satisfy the same equation, we can choose
$\overline {\boldsymbol{c}}_{\alpha ,2}^{a+} = \overline {\boldsymbol{c}}_{\alpha ,2}^{a-} = \overline {\boldsymbol{c}}_{\alpha ,2}^{a}$
; further to this, we may note that (3.8) implies that
$\overline { \unicode{x1D647}}_{\,3}^{\,a}\overline { \unicode{x1D647}}_2^a = {\overline {\boldsymbol{\varGamma }}^a}^2$
, where
$\overline {\boldsymbol{\varGamma }}^a$
is a diagonal matrix with entries given by the elements of
$\{\overline {\gamma }_{\,\alpha}^{\,a}\}_{\alpha = 0}^{\alpha _{\textrm{max}}}$
. Defining a matrix
$\overline { \unicode{x1D63E}}^a$
with each column being a different
$\overline {\boldsymbol{c}}_{\alpha ,2}^a$
(arranged so as to match the ordering of
$\overline {\boldsymbol{\varGamma }}^a$
), we then see that
\begin{equation} \overline { \unicode{x1D647}}_{\,3}^{\,a}\overline { \unicode{x1D647}}_2^a\overline { \unicode{x1D63E}}^a = \overline { \unicode{x1D63E}}^a{\overline {\boldsymbol{\varGamma }}^a}^2 = \overline { \unicode{x1D63E}}^a\overline { \unicode{x1D647}}_{\,3}^{\,a}\overline { \unicode{x1D647}}_2^a. \end{equation}
Since
$\overline { \unicode{x1D63E}}^a$
commutes with a diagonal matrix with unique entries, it too must be diagonal, so having not yet normalised the eigenvectors, we can now simply choose
$\overline { \unicode{x1D63E}}^a = \unicode{x1D644}$
. The admittance is then simple to construct from substituting it into the lower block of (3.6), getting
$\overline { \unicode{x1D647}}_{\,3}^{\,a}\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,a\pm } = \pm \overline {\boldsymbol{\varGamma }}^a$
. In both two and three dimensions,
$\overline { \unicode{x1D647}}_{\,3}^{\,a} = {\textrm{i}}\kern0.5pt a\omega \unicode{x1D644}$
, so we have
\begin{equation} \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,a\pm } = \pm \frac {1}{{\textrm{i}}\kern0.5pt a\omega }\overline {\boldsymbol{\varGamma }}^a =: \pm \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a}. \end{equation}
For clarity, the admittances in two dimensions have the following explicit expression:
\begin{equation} \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}_{\alpha \beta }^{a\pm } = \pm {\textrm{i}}\kern0.5pt\delta _{\alpha \beta }\exp\! \left \{-\frac {{\textrm{i}}\kern0.5pt \pi}{4}\left [{\textrm{sgn}}\left (1 - \frac {\lambda _\alpha ^2}{a^2\omega ^2X^2}\right ) + 1\right ]\right \}\sqrt {\left |1 - \frac {\lambda _\alpha ^2}{a^2\omega ^2X^2}\right |}, \end{equation}
and the three-dimensional case is identical, only with an
$R$
in place of every
$X$
. We can then easily write down the definition of the straight-duct cutoff frequency for each mode in two and three dimensions
\begin{equation} \overline {\omega }_\alpha ^a(s) = \begin{cases} \dfrac {\lambda _\alpha }{aX(s)} \quad &2\mbox{-}\textrm{D}, \\[9pt] \dfrac {\lambda _\alpha }{aR(s)} \quad &3\mbox{-}\textrm{D}, \\ \end{cases} \end{equation}
and note that
$\overline {\omega }_\alpha ^a$
increases with the spatial modenumber, while having the opposite relationship with the temporal modenumber.
In order to calculate the nonlinear characteristic admittances, we substitute the linear characteristic admittances into (3.2b ), and look for fixed points
\begin{align} 0 = &-\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a}\kern1.5pt \overline { \kern-1.5pt \mathcal{N}}_7^{\,\,ab}\langle \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a - b}, \unicode{x1D644}\,\rangle + \kern1.5pt \overline { \kern-1.5pt \mathcal{N}}_1^{\,\,ab}\langle \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a - b},\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,b}\rangle + \kern1.5pt \overline { \kern-1.5pt \mathcal{N}}_6^{\,\,ab} \nonumber \\&\mp \overline {\mathcal{Y}}^{\,ab\pm }\Big [\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a}\overline { \unicode{x1D647}}_{\,3}^{\,a} + \langle \overline { \unicode{x1D647}}_{\,3}^{\,a - b}\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a - b}, \unicode{x1D644}\,\rangle + \langle \unicode{x1D644},\overline { \unicode{x1D647}}_{\,3}^{\,b}\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,b}\rangle \Big ]. \end{align}
If we note that the square-bracketed quantity is actually equal to
$\overline {\boldsymbol{\varGamma }}^a + \langle \overline {\boldsymbol{\varGamma }}^{a - b}, \unicode{x1D644}\,\rangle + \langle \unicode{x1D644},\overline {\boldsymbol{\varGamma }}^b\rangle$
, a sum of three diagonal matrices each acting on a different component, then we may simply divide through by this, pointwise, for each entry of
$\overline {\mathcal{Y}}^{\,ab\pm }$
, getting
\begin{equation} \overline {\mathcal{Y}}_{\,\alpha \beta \gamma }^{\,ab\pm } = \pm \frac {\big (\kern1.5pt \overline { \kern-1.5pt \mathcal{N}}_1^{\,\,ab}\langle \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a - b},\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,b}\rangle + \kern1.5pt \overline { \kern-1.5pt \mathcal{N}}_6^{\,\,ab} - \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a}\kern1.5pt \overline { \kern-1.5pt \mathcal{N}}_{\,\,7}^{\,\,ab}\langle \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a - b}, \unicode{x1D644}\,\rangle \big )_{\alpha \beta \gamma }}{\overline {\gamma }_{\,\alpha}^{\,a} + \overline {\gamma }_{\,\beta}^{\,a - b} + \overline {\gamma }_{\,\gamma}^{\,b}} =: \pm \overline {\mathcal{Y}}_{\,\alpha \beta \gamma }^{\,ab}. \end{equation}
Obviously we will run into problems if the denominator of this turns out to be zero: this would constitute a resonant triad (see e.g. Protas, Kang & Bustamante Reference Protas, Kang and Bustamante2024). Empirically, we do not seem to run into this problem in the straight-duct case, but for the more complicated characteristic admittances that follow there seems to be no such guarantee. We avoid the problem of singular entries in the characteristic admittance by excising them from the matrix when a resonant triad is found, but proper work on the avoidance/physical significance of them should at some point be conducted.
3.2.2. Other characteristic admittances
If instead we were to consider the characteristic admittance of an infinitely long duct with constant curvature, the eigenvalues would remain on the real/imaginary axes, but the eigenvector matrices would cease to be diagonal, resulting in a slightly more complicated expression for the corresponding characteristic admittance
$\breve { \kern-3pt \unicode{x1D654}}^{\,a}$
. Including torsion in addition would cause the eigenvalues to bifurcate into the complex plane: as a result the torsional-duct characteristic admittances
$\kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{\,\,a\pm}$
would no longer satisfy
$\kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{\,\,a+} = -\kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{\,\,a-}$
. For completeness, expressions for these quantities are derived explicitly in Appendix C, although here we will only make use of the straight-duct admittance above.
3.3. Splitting operators
By using the characteristic admittances defined above, we can split waves into forward- and backward-going parts
\begin{align} &\boldsymbol{u}^a = \boldsymbol{u}^{a+} + \boldsymbol{u}^{a-}, \quad \boldsymbol{p}^{a} = \boldsymbol{p}^{a+} + \boldsymbol{p}^{a-}. \end{align}
We may perform this decomposition from the total pressure by defining linear and nonlinear splitting operators
$ \unicode{x1D64E}^{\,a\pm }$
and
$\mathcal{S}^{\,ab\pm }$
\begin{equation} \boldsymbol{p}^{a\pm } = \unicode{x1D64E}^{\,a\pm }\boldsymbol{p}^a + \sum _b\mathcal{S}^{\,ab\pm }\langle \boldsymbol{p}^{a - b},\boldsymbol{p}^b\rangle ; \end{equation}
a similar procedure could be defined for the velocity. Assuming we have picked a characteristic admittance, we now proceed to calculate these splitting operators in terms of that characteristic admittance. Restricting ourselves to the linear case for a moment, we have
\begin{equation} \boldsymbol{p}^{a\pm } = \boldsymbol{p}^a - \boldsymbol{p}^{a\mp } = \boldsymbol{p}^a - \unicode{x1D655}^{\,a\mp}(\boldsymbol{u}^a - \boldsymbol{u}^{a\pm }) = ( \unicode{x1D644} - \unicode{x1D655}^{\,a\mp} \unicode{x1D654}^{\,a})\boldsymbol{p}^a + \unicode{x1D655}^{\,a\mp} \unicode{x1D654}^{\,a\pm}\boldsymbol{p}^{a\pm }, \end{equation}
so that
\begin{equation} \boldsymbol{p}^{a\pm } = ( \unicode{x1D644} - \unicode{x1D655}^{\,a\mp} \unicode{x1D654}^{\,a\pm})^{-1}( \unicode{x1D644} - \unicode{x1D655}^{\,a\mp} \unicode{x1D654}^{\,a})\boldsymbol{p}^a = \unicode{x1D64E}^{\,a\pm }\boldsymbol{p}^a, \end{equation}
which may equivalently be written more compactly as
\begin{equation} \unicode{x1D64E}^{\,a\pm } = ( \unicode{x1D654}^{\,a\pm} - \unicode{x1D654}^{\,a\mp})^{-1}( \unicode{x1D654}^{\,a} - \unicode{x1D654}^{\,a\mp}). \end{equation}
If the above equations are expanded to second order in terms of
$ \unicode{x1D64E}^{\,a\pm }$
, we find
$\mathcal{S}^{\,ab\pm }$
to be
\begin{equation} \mathcal{S}^{\,ab\pm } = ( \unicode{x1D654}^{\,a\pm} - \unicode{x1D654}^{\,a\mp})^{-1}\big (\mathcal{Y}^{\,ab} - \mathcal{Y}^{\,ab+}\langle \unicode{x1D64E}^{\,(a - b)+}, \unicode{x1D64E}^{\,b+}\rangle - \mathcal{Y}^{\,ab-}\langle \unicode{x1D64E}^{\,(a - b)-}, \unicode{x1D64E}^{\,b-}\rangle \big ). \end{equation}
The decomposition (3.16) may worry some readers, in that our model is weakly nonlinear but here we assume a linear decomposition, i.e. that
$\boldsymbol{p}^{a} = \boldsymbol{p}^{a+} + \boldsymbol{p}^{a-}$
. This also worried the authors. A weakly nonlinear correction term to this decomposition may be added,
$\boldsymbol{p}^{a} = \boldsymbol{p}^{a+} + \boldsymbol{p}^{a-} + \sum _b\mathcal{C}^{ab}\langle \boldsymbol{p}^{a-b},\boldsymbol{p}^{b}\rangle$
, and a similar term to the
$\boldsymbol{u}^a$
decomposition. Simple algebra then shows that any arbitrary choice of these correction terms is mathematically consistent, and so here we take them to be zero; indeed, taking the cross-term to be non-zero would break the interpretation of a decomposition into upstream and downstream modes. Note that weak nonlinearity is still accounted for in the decomposition through the
$\mathcal{S}^{\,ab\pm }$
terms in (3.17). In fact, this may not be surprising when we consider that the acoustic disturbance is a sum over modes,
$p' = \sum _{a,\alpha }P^a_\alpha \psi _\alpha {\textrm{e}}^{-{\textrm{i}}\kern0.5pt a\omega t}$
, which is necessarily a linear decomposition at each axial duct location, and that the separation into downstream- and upstream-propagating sound just involves rewriting this modal expansion as two sums.
3.4. Pressure boundary condition
When we come to solve these equations, we prescribe the admittances at the duct’s outlet, amounting to a radiation condition, before solving, in order, (3.2a
) and (3.2b
). After this, we know everything about the global radiative properties of the duct, and correspondingly we have only a first-order ODE (3.3) to solve for the pressure. For musical instruments, the pressure perturbation is prescribed by the musician, who disturbs the internal air column of the instrument at the inlet, meaning that the pressure is specified at the entrance (taken here to be at
$s = 0$
) and solved forwards from there.
One matter that is not entirely clear from a mathematical point of view is that of the musician’s interaction with reflecting waves. Should the musician’s playing amount to a prescription of the forward-going pressure only, meaning that waves reflected back impact without consequence on the source? Or should it rather be a prescription of the total pressure, meaning that the musician is adapting to the impact of the reflections they themselves generated? In brass instruments, it is thought that a feedback loop is formed between backward-going waves and the player’s lip response to them: however, this complication is avoided here for simplicity, and we choose only to specify the forward-going pressure at the inlet,
$p^+(\boldsymbol{x},t)|_{s = 0}$
, with corresponding coefficient
$\boldsymbol{p}^{a+}(0)$
, and convert this (for the purpose of solving the equation) to the total pressure with the inverse splitting operator
\begin{equation} \boldsymbol{p}^a(0) = ( \unicode{x1D64E}^{\,a+})^{-1}\Big [\boldsymbol{p}^{a+}(0) - \sum _b\mathcal{S}^{\,ab+}\big \langle ( \unicode{x1D64E}^{\,(a - b)+})^{-1}\boldsymbol{p}^{(a - b)+}(0),( \unicode{x1D64E}^{\,b+})^{-1}\boldsymbol{p}^{b+}(0)\big \rangle \Big ]. \end{equation}
The simplest boundary condition for the pressure is that of a sinusoidal piston source at
$s = 0$
, given by
$p^+(\boldsymbol{x},t)|_{s = 0} = M\sin (\omega t)$
(where we recall
$M$
is the perturbation Mach number, introduced in § 2.1), and projected onto the coordinate basis as
\begin{equation} P_\alpha ^{a+}(0) = \frac {M\sqrt {A_{\textit{cs}}(0)}{\textrm{sgn}}(a)\delta ^{|a|,1}\delta _{\alpha 0}}{2{{i}}}, \end{equation}
where
$A_{{cs}}(0) = X(0)$
in two dimensions and
$ \pi R(0)^2$
in three. We can also consider non-plane modes. In two dimensions, for a variable-width straight duct (with width variation symmetric about the centreline), symmetric and antisymmetric modes are uncoupled, so there is a family of pressure distributions inaccessible to the plane-wave inlet condition. The ‘fundamental’ member of this family is the first antisymmetric mode, so we can consider a variant on condition (3.23) where the
$\delta _{\alpha 0}$
is replaced with a
$\delta _{\alpha 1}$
.
This effect is amplified for any variable-width straight duct in three dimensions, where modes of different azimuthal wavenumber are uncoupled, resulting in uncountably many mutually inaccessible families of pressure distributions. The ‘fundamental’ of each therefore represents another possible inlet condition, meaning we replace the
$\delta _{\alpha 0}$
with a
$\delta _{p_\alpha n}\delta _{q_\alpha 1}$
for any
$n$
of our choice.
3.5. Numerical method
3.5.1. Truncation
When solving this infinite set of coupled ODEs numerically, we must truncate both the number of spatial and temporal modes used. If we pick maximum values of each, we then have spatial modes ranging from
$0$
to
$\alpha _{\textrm{max}}$
and temporal modes ranging from
$-a_{\textrm{max}}$
to
$a_{\textrm{max}}$
. This would mean that we have
$(\alpha _{\textrm{max}}+1)(2a_{\textrm{max}}+1)$
equations to solve for
$P_\alpha ^a$
,
$(\alpha _{\textrm{max}}+1)^2(2a_{\textrm{max}}+1)$
for
$ \unicode{x1D654}_{\,\alpha \beta }^{\,a}$
and
$(\alpha _{\textrm{max}}+1)^3(2a_{\textrm{max}}+1)^2$
for
$\mathcal{Y}_{\,\alpha \beta \gamma }^{\,ab}$
. However, by our previous assumptions we are neglecting zero modes, and all of our acoustic quantities being real implies that
$P^{-a}_\alpha = P^{a*}_\alpha$
and
$U^{-a}_\alpha = U^{a*}_\alpha$
, resulting in
$(\alpha _{\textrm{max}}+1)a_{\textrm{max}}$
independent equations for
$P_\alpha ^a$
and
$(\alpha _{\textrm{max}}+1)^2a_{\textrm{max}}$
for
$ \unicode{x1D654}_{\,\alpha \beta }^{\,a}$
.
A slightly messier calculation at this point puts the number of independent entries of
$\mathcal{Y}_{\,\alpha \beta \gamma }^{\,ab}$
at
$(\alpha _{\textrm{max}}+1)^3(2a_{\textrm{max}}-1)a_{\textrm{max}}$
. More entries may yet be discarded, however, since discrete convolutions behave strangely upon truncation; if
$a$
and
$b$
range from
$-a_{\textrm{max}}$
to
$a_{\textrm{max}}$
then
$a {-} b$
ranges from
$-2a_{\textrm{max}}$
to
$2a_{\textrm{max}}$
, so as a consequence, at each time step, for every triple
$(\alpha ,\beta ,\gamma )$
, there is an
$ab$
matrix
$\mathcal{Y}_{\,\alpha \beta \gamma }^{\,ab}$
of which the only calculable entries form a banded matrix (with bandwidth
$a_{\textrm{max}}$
), resulting in
$a_{\textrm{max}}(a_{\textrm{max}}+1)$
lost entries. Having exploited the conjugacy property already, this means a subtraction of
$a_{\textrm{max}}(a_{\textrm{max}}+1)/2$
entries from the total, resulting in
$3(\alpha _{\textrm{max}}+1)^3a_{\textrm{max}}(a_{\textrm{max}}-1)/2$
independent entries for the nonlinear admittance.
3.5.2. Numerical solver
We use a fourth-order Runge–Kutta method to integrate the admittance and pressure equations in
$s$
. Unless otherwise specified, we use the adaptive-step solver ode45 in Matlab, although for one case in § 4.4 we use a fixed-step solver, as that proves more accurate in the neighbourhood of a cusp in the duct diameter.
3.5.3. Numerical viscosity
One of the pitfalls of truncation is the pooling of energy at higher modes. This can be countered by employing a numerical viscosity, which we do by subtracting a new term
$E^a$
from the right-hand-side of the pressure (3.3). Here, we take this to be
\begin{equation} E^a = -\nu _0\frac {|a|\omega \beta_0 M}{1 + M \beta_0\omega s}\log\! \left (1 - \frac {|a| - 1}{a_{\textrm{max}}}\right )P_0^a \sim \nu _0\frac {|a|^2\omega \beta_0 M}{a_{\textrm{max}}}P_0^a, \end{equation}
where
$\nu _0$
is a positive scalar determining how great an effect we want the viscosity to have. The second approximation here is in the limit of small
$s$
relative to the shock-formation distance
$1/M \beta_0\omega$
, and with
$1 \ll |a| \ll a_{\textrm{max}}$
(where energy pooling starts to occur). This choice of numerical viscosity is informed by the truncation error for a sawtooth plane wave in a one-dimensional duct, and is discussed further in § 4.1 below.
The physical molecular viscosity could also instead be considered as a starting point, and comparing the result of one-dimensional numerical sawtooth stabilisation with its physical counterpart, we find that the first-order correction term agrees with the
$a^2$
factor from the second approximation, but also that the combined shear and bulk viscosities
$4\mu /3 + \zeta$
must be approximately equal to
$2 \beta_0 M/a_{\textrm{max}}\omega$
. This calculation may be found in Appendix D.
3.5.4. Alternative numerical damping
Numerical instabilities can also arise in linear problems. If a node occurs in the pressure (scenarios involving this are discussed in § 4.5), this causes a singularity in the admittance. Such singularities may be avoided (or ‘dampened’) by the addition of a small imaginary part to the frequency
$\omega$
, which corresponds to a sound source slowly exponentially growing in time, and therefore a wave pattern that slowly exponentially decays in space away from the source. This technique is borrowed from stability analysis (e.g. Briggs Reference Briggs1964; Bers Reference Bers1983), where adding an imaginary part to the frequency
$\omega$
is used to guarantee a causal solution.
3.6. Modal resolution
The following test cases require different numbers of modes depending on what physics is being considered. For linear acoustics, spatial modes are the only limiting factor: when comparing with pre-existing linear work, we use the previous resolution as a benchmark (e.g.
$\alpha _{\textrm{max}} = 30$
for figure 4), but when comparing with results, we pick a new standard (taking
$\alpha _{\textrm{max}} = 50$
for figures 7, 8, 9, 10). In all of these cases the absence of nonlinearity allows us to set
$a_{\textrm{max}} = 1$
, and unoptimised Matlab code run on a standard laptop finishes in the order of seconds.
Where nonlinearity comes in, memory requirements are likely to exceed that of a standard laptop, although all results here have been computed on a single desktop with
$128\,\textrm{GiB}$
of memory. For most of these computations (i.e. Figures 3 and 6–16) spatial resolution is necessary, so we take
$a_{\textrm{max}} = 10$
and
$\alpha _{\textrm{max}} = 10$
, and computations take of the order of tens of minutes. For figure 2, we consider only plane waves in a straight duct, so we can take
$\alpha _{\textrm{max}} = 0$
and
$a_{\textrm{max}} = 100$
; while this does not exceed the memory limit, the calculation now takes of the order of half a day to compute using a single core.
The numerical procedure converges as the number of temporal and spatial modes is increased, and the truncations above are sufficient for accurate results in all of the cases presented here. A demonstration of the numerical convergence with the number of modes is provided for a constant-curvature bend in three dimensions in § 4.3 below.
4. Results
We now consider applying the procedure above to a number of test cases. The numerical Matlab code to compute these cases is provided in the supplementary material.
4.1. A straight duct of constant diameter
The simplest test geometry in either two or three dimensions is that of a straight duct with constant diameter. For a plane piston source, the equations are simplified greatly, both by the absence of curvature and width-variation terms, and also by the collapsing of the spatial-mode vectors onto a single scalar component. The admittance equations are then solved trivially, since the boundary condition is a constant solution, giving
\begin{align} \unicode{x1D654}_{00}^a(s) = 1, \quad \mathcal{Y}_{000}^{ab}(s) = -\frac {1}{2\sqrt {A_{{cs}}}}\left ( \beta_0 + \frac {a - 2b}{a}\right )\quad \forall s. \end{align}
This holds in both two and three dimensions, and we write
$A_{{cs}}$
for the area of the duct cross-section so that we can work in generality. Plugging these into the pressure equation, we get (again in both cases, once we substitute the value of
$\mathcal{I}_{000}$
)
\begin{equation} \frac {{\rm d} P_0^a}{{\rm d} s} = {\textrm{i}}\kern0.5pt a \omega P_0^a + \frac {{\textrm{i}}\kern0.5pt\omega }{2\sqrt {A_{{cs}}}}\sum _b\big [-a \beta_0 + (a - 2b)\big ]P_0^{a - b}P_0^b, \end{equation}
but the second term in the square bracket here vanishes upon summation. Since the zero mode
$\psi _0$
does not depend on
$s$
for a straight duct, we can multiply through by it, getting an equation for Fourier coefficients
\begin{equation} \frac {{\rm d} P^a}{{\rm d} s} = {\textrm{i}}\kern0.5pt a\omega P^a - \frac {{\textrm{i}}\kern0.5pt a\omega \beta_0}{2}\sum _bP^{a - b}P^b. \end{equation}
This equation was solved up to shock formation by Fubini (Reference Fubini1935) (for any boundary conditions both periodic and odd in t), in the form of a sine series for the pressure
\begin{align} &\frac {p}{M} = \sum _aB_a\sin [a\omega (t - s)];&& B_a = \frac {2}{a\sigma }{{J}}_a(a\sigma ), \end{align}
where
$\sigma = M \beta_0\omega s$
is the arclength normalised by the shock-formation distance. Fay (Reference Fay1931) found a post-shock-formation solution in the form of a sawtooth wave
\begin{equation} B_a = \frac {2}{a(1 + \sigma )}, \end{equation}
valid for
$\sigma \gtrapprox 3$
. Blackstock (Reference Blackstock1966) matched these two solutions with the following two terms:
\begin{equation} B_a = \frac {2}{a \pi}P_{\textit{sh}} + \frac {2}{a \pi\sigma }\int _{\varPhi _{\textit{sh}}}^ \pi\cos a(\varPhi - \sigma \sin \varPhi )\,{\rm d}\varPhi , \end{equation}
where the shock amplitude
$P_{\textit{sh}}$
and phase
$\varPhi _{\textit{sh}}$
satisfy
\begin{align} &P_{\textit{sh}} = \sin (\sigma P_{\textit{sh}}),&& \varPhi _{\textit{sh}} = \sigma \sin \varPhi _{\textit{sh}}. \end{align}
For
$\sigma \lt 1$
, the only solution for
$(P_{\textit{sh}},\varPhi _{\textit{sh}})$
is
$(0,0)$
, so the first term vanishes and the second becomes the Fubini solution. As
$\sigma$
grows greater than 1, two more solutions appear either side of 0; the relevant solutions here are the positive ones. For
$\sigma \gtrapprox 5 \pi/2$
more solutions appear; we remain interested from this point in the minimal positive solution for each of
$P_{\textit{sh}}$
and
$\varPhi _{\textit{sh}}$
. Since
$\varPhi _{\textit{sh}}$
converges to
$ \pi$
from below, the second term tends to 0, while the first term behaves like the Fay solution.
The numerical viscosity used here, described in § 3.5.3 above, is justified by considering the truncation error of a sawtooth wave when substituted into (4.3). The sawtooth wave as derived in Fay (Reference Fay1931) has Fourier coefficients given by
\begin{equation} P^a = \frac {{{i}} \textit{Me}^{{\textrm{i}} a\omega s}}{a(1 + \sigma )} , \end{equation}
so the terms missing if we truncate the nonlinear term in (4.3) to a finite sum are then
\begin{equation} \begin{aligned} E^a &= \frac {{\textrm{i}} a\omega \beta_0 \sqrt {A_{{cs}}}}{2}\sum _{\substack {b = -\infty ,\\b \neq 0,a}}^\infty P^{a - b}P^b - \frac {{\textrm{i}} a\omega \beta_0 \sqrt {A_{{cs}}}}{2}\sum _{\substack {b = a - {\textrm{sgn}}(a)a_{\textrm{max}},\\b \neq 0,a}}^{{\textrm{sgn}}(a)a_{\textrm{max}}}P^{a - b}P^b \\ &= -\frac {{\textrm{i}}\kern0.5pt\omega \beta_0 \sqrt {A_{{cs}}}M^2\textrm{e}^{{\textrm{i}}\kern0.5pt a\omega s}}{(1 + \sigma )^2}\left (\sum _{\substack {b = -\infty ,\\b \neq 0,a}}^\infty - \sum _{\substack {b = a - {\textrm{sgn}}(a)a_{\textrm{max}},\\b \neq 0,a}}^{{\textrm{sgn}}(a)a_{\textrm{max}}}\right )\frac {a}{2(a - b)b} ,\end{aligned} \end{equation}
where the limits on the second sum ensure that for both positive and negative
$b$
,
$|a - b|$
may never exceed
$a_{\textrm{max}}$
. By carefully splitting the summands into partial fractions, which make the sums easier to compute, we gain a form which has integral bounds that we use to calculate the numerical viscosity. This is done in detail in Appendix D, and results (with the inclusion of the viscous scale factor
$\nu _0$
) in (3.24).
Comparison of numerically calculated mode amplitudes (dashed) with those predicted by the Blackstock solution (solid), using the code in two dimensions. (a) Shows modes 1 through 7, while (b) shows 10 through 70. Truncation was taken at
$\alpha _{\textrm{max}} = 0$
and
$a_{\textrm{max}} = 100$
.
Pressure (normalised to the source amplitude) in a bend identical to the one used in McTavish & Brambley (Reference McTavish and Brambley2019), with a plane piston source of frequency
$3/X$
, for (a) linear, (b)
$M$
= 0.05, (c) 0.10 and (d) 0.15. Truncation was taken at
$\alpha _{\textrm{max}} = a_{\textrm{max}} = 10$
. An animated version of this figure is available as Movie 1 in the supplementary material.
Pressure (normalised to the source amplitude) in a bend identical to the one used in Félix & Pagneux (Reference Félix and Pagneux2002), with a plane piston source of frequency
$2.4/R$
, plotted both through the midplane and across the duct outlet in the linear regime. Truncation was taken at
$\alpha _{\textrm{max}} = 30$
and
$a_{\textrm{max}} = 1$
. An animated version of this figure is available as Movie 2 in the supplementary material.
Numerical convergence in a three-dimensional planar bend for three definitions of the error, detailed in (4.10). Here,
$a_{\textrm{max}}=1$
, and so the number of modes is
$\alpha _{\textrm{max}}+1$
.
Pressure (normalised to the source amplitude) in a bend identical to the one used in Félix & Pagneux (Reference Félix and Pagneux2002), with a plane piston source of frequency
$2.4/R$
, for (a) linear, (b)
$M$
= 0.02, (c) 0.05 and (d) 0.10. Truncation was taken at
$\alpha _{\textrm{max}} = a_{\textrm{max}} = 10$
. An animated version of this figure is available as Movie 3 in the supplementary material.
We can compare our numerics with this. This can be done in both two and three dimensions; figure 2 compares various modes in the two-dimensional case, for a simulation with 1 spatial mode (since no spatial coupling is induced for this geometry anyway) and 100 temporal modes. As one would expect, good matches are achieved for the lower (dominant) modes, while for higher modes the numerical solutions undershoot at first and over-compensate later on. We have achieved matches of equal quality with the code in three dimensions.
Eliminating spatial variation in favour of higher temporal modes is a useful special case because it allows for the calibration of the numerical viscosity. Figure 2 was created with
$\nu _0 = 1$
, because higher scale factors than this place far too much damping on the higher modes, while lower factors allow instabilities to build up much more quickly. In light of this, all other nonlinear calculations take
$\nu _0 = 1$
as well.
4.2. Constant curvature in two dimensions
In two dimensions, McTavish & Brambley (Reference McTavish and Brambley2019) considered linear and nonlinear propagation around a constant-width bend, first computed in the linear case by Félix & Pagneux (Reference Félix and Pagneux2001). This provides a test case for the combination of curvature and nonlinearity. Figure 3 reproduces figure 7 of McTavish & Brambley (Reference McTavish and Brambley2019), with all of the same parameters, i.e. (for constant duct width
$X$
) a bend of curvature
$8/5X$
located at a distance
$2X$
downstream of a plane piston source of frequency
$3/X$
, with truncation taken at 10 spatial modes and 10 temporal modes, plotted in the linear regime and nonlinearly for
$M = 0.05$
, 0.10 and 0.15. Good agreement is observed.
4.3. Constant curvature in three dimensions
In three dimensions, no published work exists that can be used to test the combined effects of curvature and nonlinearity, but in the linear case Félix & Pagneux (Reference Félix and Pagneux2002) provide an example that may be tested against. Here, for constant duct radius
$R$
, a plane piston source of frequency
$2.4/R$
is placed at the mouth of a bend of curvature
$4/5R$
, and a cross-section through the half-plane is plotted alongside the outlet’s circular cross-section. Comparing figure 4 with figure 2 of Félix & Pagneux (Reference Félix and Pagneux2002), we again see good agreement. Truncation was taken at 30 spatial modes to match their calculation, with the usual 1 temporal mode in the absence of nonlinearity.
For this geometry, we also follow Félix & Pagneux (Reference Félix and Pagneux2002) in the inclusion of a convergence study for the spatial modes. As well as their two definitions of the error involving integration along a surface, we also calculate a volumetric error, which behaves similarly. The reference solution
$p_{\textit{ref}}$
was calculated here with
$\alpha _{\textrm{max}} = 50$
and
$a_{\textrm{max}} = 1$
. Figure 5 shows the result. As with Félix & Pagneux (2002), there is uneven behaviour due to the different effects of higher azimuthal or radial resolution; nonetheless, for error metrics
$\epsilon _1$
and
$\epsilon _3$
the inclusion of more modes never causes the error to grow, and even for error metric
$\epsilon _2$
the solution is clearly still converging as more modes are included. The error definitions are
\begin{align} \epsilon _1 &= \sqrt {\frac {\int _{t = 0}^{2 \pi/\omega }\int _{s = 0}^{s_d}\int _{r = 0}^R\left [h_s\|p - p_{\textit{ref}}\|^2\right ]_{\theta = 0}^{ \pi}{\rm d}r{\rm d}s{\rm d}t}{\int _{t = 0}^{2 \pi/\omega }\int _{s = 0}^{s_d}\int _{r = 0}^R\left [h_s\|p_{\textit{ref}}\|^2\right ]_{\theta = 0}^{ \pi}{\rm d}r{\rm d}s{\rm d}t}}, \\[-12pt] \nonumber \end{align}
\begin{align} \epsilon _2 &= \sqrt {\frac {\int _{t = 0}^{2 \pi/\omega }\int _{\theta = 0}^{2 \pi}\int _{r = 0}^Rr\|p - p_{\textit{ref}}\|^2|_{s = s_d}{\rm d}r{\rm d}\theta {\rm d}t}{\int _{t = 0}^{2 \pi/\omega }\int _{\theta = 0}^{2 \pi}\int _{r = 0}^Rr\|p_{\textit{ref}}\|^2|_{s = s_d}{\rm d}r{\rm d}\theta {\rm d}t}}, \\[-9pt] \nonumber \end{align}
\begin{align} \epsilon _3 &= \sqrt {\frac {\int _{t = 0}^{2 \pi/\omega }\int _{\theta = 0}^{2 \pi}\int _{r = 0}^R\int _{s = 0}^{s_d}rh_s\|p - p_{\textit{ref}}\|^2|_{s = s_d}{\rm d}s{\rm d}r{\rm d}\theta {\rm d}t}{\int _{t = 0}^{2 \pi/\omega }\int _{\theta = 0}^{2 \pi}\int _{r = 0}^R\int _{s = 0}^{s_d}rh_s\|p_{\textit{ref}}\|^2|_{s = s_d}{\rm d}r{\rm d}s{\rm d}\theta {\rm d}t}}. \\[9pt] \nonumber \end{align}
In figure 6 the same geometry has a plane piston source of frequency
$\omega = 2.4/R$
prescribed in the total pressure at the top left, for various Mach numbers. We see good agreement both with the published linear work of Félix & Pagneux (Reference Félix and Pagneux2002), where evenly spaced peaks and troughs travel around the outside of the bend, and also with the unpublished nonlinear work of McTavish (Reference McTavish2018), where the contours are deformed by the steepening of the peaks as the Mach number is increased. To match the latter calculation, 10 spatial modes and 10 temporal modes were used.
4.4. An exponential horn
The exponential horn is a convenient geometry: if an approximation is made allowing only plane waves to propagate, an analytical solution (Webster Reference Webster1919) exists, with the growth rate of the horn acting to dampen oscillations. The geometry is best parametrised across both two and three dimensions by specifying the growth exponent of the horn’s cross-sectional area, i.e. by setting
\begin{equation} A_{{cs}}(s) = \begin{cases} A_{{cs}}(0)\exp (2ms)\quad &s\in [0,s_{{o}}],\\ A_{{cs}}(0)\exp (2ms_{{o}}) &s\gt s_{{o}}, \end{cases} \end{equation}
so that the radial growth factor in two dimensions is doubled relative to that in three. The section of duct with
$s \gt s_{{o}}$
forms the infinite straight duct providing the characteristic admittance boundary condition
$ \unicode{x1D654}^{\,a} = \unicode{x1D654}^{\,a+}$
at the outlet
$s=s_0$
. Since we are only interested in plane waves, the (scalar) linear admittance equation (in the variable radius region) becomes
\begin{equation} \frac {{\rm d} Y^a}{{\rm d} s} = -{\textrm{i}}\kern0.5pt a\omega (Y^a)^2 - 2mY^a + {\textrm{i}}\kern0.5pt a\omega , \end{equation}
with the boundary condition at
$s = s_{{o}}$
being
$Y^a = 1$
(which also holds for all
$s \gt s_{{o}}$
). Solving this ODE, we find that
\begin{equation} Y^a = \frac {n_a\cos n_a(s_{{o}} - s) - ({\textrm{i}}\kern0.5pt a\omega - m)\sin n_a(s_{{o}} - s)}{n_a\cos n_a(s_{{o}} - s) - ({\textrm{i}}\kern0.5pt a\omega + m)\sin n_a(s_{{o}} - s)}, \end{equation}
for
$s \leqslant s_{{o}}$
. Here
$n_a$
is a quantity satisfying
$n_a^2 = a^2\omega ^2 - m^2$
, which is real in the case of cut-on modes and imaginary otherwise. Having found
$Y^a$
, it may now be plugged into the equation for the pressure, which in the plane-wave case reduces to
\begin{equation} \frac {{\rm d}P_0^a}{{\rm d}s} = \begin{cases} (m + {\textrm{i}}\kern0.5pt a\omega Y^a)P_0^a \quad &s \leqslant s_{{o}}, \\ {\textrm{i}}\kern0.5pt a\omega P_0^a &s \gt s_{{o}}, \end{cases} \end{equation}
with our initial condition being
\begin{equation} P_0^a = \frac {M\sqrt {A_{{cs}}(0)}}{2{{i}}}\delta ^{|a|1}{\textrm{sgn}}(a), \end{equation}
or alternatively
\begin{equation} P_0^a = \frac {M\sqrt {A_{{cs}}(0)}}{2{{i}}}\delta ^{|a|1}{\textrm{sgn}}(a)\frac {n_a\cos n_as_{{o}} - ({\textrm{i}}\kern0.5pt a\omega + m)\sin n_as_{{o}}}{n_a\cos n_as_{{o}} - {\textrm{i}}\kern0.5pt a\omega \sin n_as_{{o}}} \end{equation}
if we wish to specify only the forward-going pressure. Solving this equation, with condition (4.15), the pressure is
\begin{equation} P_0^a(s) = \frac {M{\textrm{sgn}}(a)\delta ^{|a|,1}}{2{{i}}}\sqrt {A_{{cs}}(0)}\frac {n_a\cos n_a(s_{{o}} - s) - ({\textrm{i}}\kern0.5pt a\omega + m)\sin n_a(s_{{o}} - s)}{n_a\cos n_as_{{o}} - ({\textrm{i}}\kern0.5pt a\omega + m)\sin n_as_{{o}}}, \end{equation}
and if we specify the forward-going pressure instead,
$m$
disappears from the bracketed term in the denominator.
To test our code on this geometry, we follow McTavish & Brambley (Reference McTavish and Brambley2019) in considering a two-dimensional horn of inlet width
$X_{\textrm{i}}$
, length
$4.5X_{\textrm{i}}$
and width increase ratio of 16, and a plane piston source of frequency of
$0.95\overline {\omega }_1^1(s_{{o}}) = 0.95 \pi/X(s_{{o}})$
for the inlet condition on the total pressure (where we recall the definition of the cutoff frequency from (3.13)). With 50 spatial modes and 1 temporal mode (since with no nonlinearity there is no temporal coupling) a good match is observed with their results in the linear case (figure 7).
Further to this, a more direct comparison with Webster’s analytical approximation is shown in figure 8, where the root-mean-square (r.m.s.) pressure at the centreline is plotted (in contrast to McTavish & Brambley (Reference McTavish and Brambley2019), the analytical solution is calculated for the whole domain, including the straight-duct section post-outlet). We achieve a very good match by restricting our calculation to plane waves only: this is achieved with a fixed (and very small) stepsize method, since the discontinuity in
$X'$
at the outlet is severe enough to cause errors in a variable stepsize solver like ode45. In contrast, the 50-mode calculation disagrees very obviously with the Webster approximation, showing that the Webster approximation will necessarily fail to capture the spatial coupling effects of a duct bell.
Pressure (normalised to the source amplitude) in an exponential horn identical to the one used in McTavish & Brambley (Reference McTavish and Brambley2019), with a plane piston source of frequency
$0.95\overline {\omega }_1^1(s_{{o}}) = 0.95 \pi/X(s_{{o}})$
, for (a) linear, (b) linear, quarter of a cycle later. Truncation was taken at
$\alpha _{\textrm{max}} = 50$
and
$a_{\textrm{max}} = 1$
. An animated version of this figure is available as Movie 4 in the supplementary material.
The r.m.s. pressure along the centreline of the exponential horn for two different modal resolutions, compared with the plane-wave approximation (Webster Reference Webster1919). Mode coupling induced by the variation in duct width causes the plane-wave approximation to be an inaccurate one for this geometry.
4.5. An inverse exponential horn in two dimensions
If we consider an ‘inverse’ exponential horn, that is, one that narrows rather than widening, we introduce the possibility of nodes. This is because waves are now reflecting to a much greater degree, so cancellations can occur. This effect is most readily observed with antisymmetric waves (and in fact cannot occur in a straight duct when considering only plane waves). Since antisymmetric waves are uncoupled from symmetric ones in a straight duct, there are grounds to consider the derivation of an analogous Webster-like solution to the first antisymmetric mode. We will solve only the two-dimensional problem here, since the three-dimensional problem is complicated by a greater number of ways to break azimuthal symmetry. As with the Webster solution, this discussion also covers only linear acoustics in this geometry.
When the Webster solution is derived, we consider a single-mode truncation of (2.50) for straight ducts of variable width (whose centreline is always in the middle of the duct, so
$X_+ = - X_- = X/2$
) in the linear regime. As such, we have no
$\mathcal{N}^{\,ab}$
, and
$ \unicode{x1D647}^a$
(which for
$a = 1$
we call
$ \unicode{x1D647}$
) reads
\begin{equation} \unicode{x1D647} = \begin{pmatrix} -\dfrac {X'}{2X}\left ( \unicode{x1D652} - \kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}\right ) &{\textrm{i}}\kern0.5pt\omega \left ( \unicode{x1D644} - \dfrac {\boldsymbol{\varLambda }^2}{\omega ^2X^2}\right ) \\[8pt] {\textrm{i}}\kern0.5pt\omega \unicode{x1D644} &\dfrac {X'}{2X}\left ( \unicode{x1D652}^T - \kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}^T\right ) \end{pmatrix}. \end{equation}
The non-coupling between symmetric and antisymmetric waves is encoded here by the fact that the matrix
$ \unicode{x1D652} - \kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}$
has non-zero entries only where the column index and row index are either both even (antisymmetric) or both odd (symmetric). In particular, this results in the following 2-mode truncation of (2.50):
\begin{equation} \dfrac {{\rm d}}{{\rm d} s}\begin{pmatrix} u_0 \\ u_1 \\ p_0 \\ p_1 \end{pmatrix} = \begin{pmatrix} -\dfrac {X'}{2X}\begin{pmatrix} 1 &0 \\ 0 &2 \end{pmatrix} &{\textrm{i}}\kern0.5pt\omega \begin{pmatrix} 1 &0 \\ 0 &1 - \dfrac { \pi^2}{\omega ^2X^2} \end{pmatrix} \\[12pt] {\textrm{i}}\kern0.5pt\omega \begin{pmatrix} 1 &0 \\ 0 &1 \end{pmatrix} &\dfrac {X'}{2X}\begin{pmatrix} 1 &0 \\ 0 &2 \end{pmatrix} \end{pmatrix}\begin{pmatrix} u_0 \\ u_1 \\ p_0 \\ p_1 \end{pmatrix}. \end{equation}
Solving the symmetric problem
$(u_0,p_0)$
results in the Webster horn equation, while solving the antisymmetric problem
$(u_1,p_1)$
results in the following equation for
$p_1$
(from which we henceforth drop the subscript 1, as with
$u_1$
):
\begin{equation} \frac {{\rm d}^2p}{{\rm d} s^2} - \left [4m^2 - \omega ^2 + \frac { \pi^2}{X_{\textrm{i}}^2}{\textrm{e}}^{4|m|s}\right ]p = 0, \end{equation}
where for inlet width
$X_{\textrm{i}}$
we have now substituted in the two-dimensional inverse exponential horn definition
$X := X_{\textrm{i}}{\textrm{e}}^{-2|m|s}$
. If we use the substitution
$\sigma = \pi{\textrm{e}}^{2|m|s}/2|m|X_{\textrm{i}}$
, this reduces to the modified Bessel equation
\begin{equation} \sigma ^2\frac {{\rm d}^2p}{{\rm d} \sigma ^2} + \sigma \frac {{\rm d} p}{{\rm d}\sigma } - \big(\nu ^2 + \sigma ^2 \big)p = 0, \end{equation}
where
$\nu$
is defined by
$\nu ^2 = 1 - \omega ^2/4m^2$
, and takes purely real values or purely imaginary values, depending on whether
$\omega$
exceeds
$2|m|$
. This equation has two solutions,
${{I}}_\nu (\sigma )$
and
${{K}}_\nu (\sigma )$
. In Appendix E we derive the solution for the pressure in detail, finding that
${{K}}_\nu (\sigma )$
dominates for most of the duct. Nodes are then found at roots of this function, i.e. roots
$s_{{node}}$
of
${{K}}_{{\textrm{i}}\sqrt {\omega ^2/4m^2 - 1}}( \pi{\textrm{e}}^{2|m|s}/2|m|X_{\textrm{i}})$
.
We can test our code on this geometry as well. To form a good physical picture of where the admittance singularities are, we may exploit the lack of coupling between symmetric and antisymmetric modes in this geometry to eliminate the plane-wave component from the solution altogether. To do this we make use of the alternative inlet condition mentioned in § 3.4. Figure 9 shows the pressure distribution lined up with a demonstration of the admittance singularity in this geometry. This was calculated with 50 spatial modes and 1 temporal mode (since we are testing just linearity here), for a duct of length
$4.5X_{\textrm{i}}$
and decrease ratio of 4. As mentioned in § 3.5.4 we add a small imaginary part to the frequency to ‘dampen’ any singularities. The frequency is then
$(5 + 0.01{\textrm{i}})/X_{\textrm{i}}$
, resulting in a value of
$\textrm{Im}(\nu )$
high enough that the first root of
${{K}}_\nu (\sigma )$
is within the duct domain, causing a node in the pressure, observable a little beyond
$s = 0.5$
.
Pressure (top, normalised to the source amplitude) inside an inverse exponential horn, with an antisymmetric source of frequency
$(5 + 0.01{{i}})/X_{\textrm{i}}$
, together with the normalised deviation of the admittance from the characteristic admittance (bottom). Truncation was taken at
$\alpha _{\textrm{max}} = 50$
and
$a_{\textrm{max}} = 1$
. An animated version of this figure is available as Movie 5 in the supplementary material.
Figure 10 compares the results of a single-spatial-mode calculation with the 50-mode calculation (the same frequency, and therefore the same stabilising imaginary part, were used on the single-mode calculation). We find that the coupling is very minimal: while the 1-mode line matches the analytical solution exactly, the 50-mode line is a lot closer than it was in the case of the Webster plane-wave solution. We conclude from this that antisymmetric sources induce less coupling than symmetric plane-wave ones.
The r.m.s. pressure along the wall of the exponential horn for two different modal resolutions, compared with the analytical solution.
4.6. Curvature and width variation combined
We can use the inverse exponential horn to demonstrate mode coupling. In the linear case, if we send in an antisymmetric source and the duct is straight, only antisymmetric modes will be excited. This symmetry can be broken in multiple ways. If nonlinearity is included, higher-order temporal harmonics will be coupled, and this can include symmetric modes (a physical argument for this is that the nonlinear terms are quadratic in the pressure, so two odd pressures will multiply to create an even one). Alternatively symmetric modes of the same temporal order can be picked up if the originally straight duct is bent. Symmetric modes can include plane waves, which have the property of always being cut-on, so we see that for an antisymmetric source of a certain frequency, propagating waves will only ‘escape’ the duct if (a) they are allowed to steepen or (b) the duct is bent.
Figure 11 displays this phenomenon of ‘acoustic leakage’, for the same inverse exponential horn as in figure 9. An antisymmetric inlet source of frequency
$1/X_{\textrm{i}}$
is used: this ensures that all non-plane waves are cutoff at the outlet. We see plane-wave tunnelling induced both for a straight duct with Mach number 0.05, or in the linear regime with a curvature of
$0.2/X_{\textrm{i}}$
. Truncation was taken at 10 spatial modes and 10 temporal modes for each of these calculations.
Pressure (normalised to the source amplitude) for an antisymmetric source of frequency
$1/X_{\textrm{i}}$
in an inverse exponential horn for (a) linear, straight, (b)
$M=0.05$
, straight and (c) linear,
$\kappa = 0.2/X_{\textrm{i}}$
. Truncation was taken at
$\alpha _{\textrm{max}} = a_{\textrm{max}} = 10$
. An animated version of this figure is available as Movie 6 in the supplementary material.
While acoustic leakage from an inverse exponential horn does not bear particular application to brass instruments, it is an interesting object of study in its own right given the potential wider scope of duct acoustics beyond musical instruments. It is included here as example of a phenomenon that can now be analysed owing to the present framework.
4.7. Torsion
Our framework allows us to investigate the effects of torsion on both linear and nonlinear acoustics. To this end, we examine helical ducts of constant wall radius
$R$
and curvature
$\kappa$
, with various torsions
$\tau$
. For each helix, we have a duct centreline given by
\begin{equation} \boldsymbol{q}(s) = (a\cos \tilde {s},a\sin \tilde {s},b\tilde {s}), \end{equation}
where
$a$
is the helical radius,
$2 \pi b$
is the pitch of the helix, and
$\tilde {s}$
is a transformed arclength. These are given in terms of
$\kappa$
and
$\tau$
by
\begin{align} a = \frac {\kappa }{\kappa ^2 + \tau ^2}, &&b = \frac {\tau }{\kappa ^2 + \tau ^2}, &&\tilde {s} = \sqrt {\kappa ^2 + \tau ^2}s. \end{align}
The binormal is calculated as
$\boldsymbol{b} = \sqrt {\kappa ^2 + \tau ^2}(b\sin \tilde {s},-b\cos \tilde {s},a)$
. At each point along the centreline, this vector points upward with a degree of backward tilt parallel to the centreline to make it perpendicular to the tangent vector. The angle between
$\boldsymbol{b}$
and the vertical is calculated with the dot product to be
\begin{equation} \cos \theta _{\boldsymbol{b}} = \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{e}_z = a\sqrt {\kappa ^2 + \tau ^2}, \end{equation}
so if we wish to know the vertical distance from a point
$\tilde {s}_0$
to the duct wall above, we consider the hypotenuse of a right-angled triangle with adjacent side of length
$R$
and an angle of
$\theta _{\boldsymbol{b}}$
between. This distance is therefore
\begin{equation} \frac {R}{\cos \theta _{\boldsymbol{b}}} = \frac {R\sqrt {\kappa ^2 + \tau ^2}}{\kappa }. \end{equation}
We also know from the definition of the helix that the vertical distance between point
$\tilde {s}_0$
and point
$\tilde {s}_0 + 2 \pi$
(where the centreline next passes over it) is
$2 \pi b$
. Therefore, for non-self-intersection, we require that
$2 \pi b$
be greater than twice the vertical centreline-to-wall distance calculated above. This constraint yields the following condition for the helical duct:
\begin{equation} 0 \gt (\tau R)^6 + 3(\kappa R)^2(\tau R)^4 + \big (3(\kappa R)^2 - \pi^2\big )(\kappa R)^2(\tau R)^2 + (\kappa R)^6. \end{equation}
All of the helical ducts we consider have curvature given by
$\kappa R = 2/3$
. We consider three torsions:
$\tau R$
= 0.16 (just large enough to avoid self-intersection),
$\tau R = 0.2$
, and
$\tau R = 1$
. All are compared in the linear case in figure 12, while the second is held fixed for varying Mach number in figure 13. We use a frequency of
$\omega = 0.95\overline {\omega }_1^1 = 0.95 \times 1.8412/R$
and truncation for all helical calculations is taken at 10 spatial modes and 10 temporal modes.
Each of these examples shows how torsional coupling causes a planar source to become non-planar at the duct outlet. It is also clear from figure 12 that (for all of the torsions we consider) the highest peaks and lowest troughs in the pressure occur on the outside of the helix, which is consistent with what was observed both in two- and three-dimensional planar bends. Figure 13 shows that when the amplitude is high enough, the peaks and troughs then lose their even spacing as the wave steepens along the outside of the helix. The effect of the
$\theta _0(s)$
coordinate twisting due to the helical coordinate system is less visible in the static images shown in figures 12 and 13, but is more evident in the animated versions of them shown in Movies 7 and 8 in the supplementary material; the wave at the outlet is seen to spin as it propagates due to the helical nature of the duct.
Linear pressure field (normalised by the Mach number) on the walls of a helical duct for the three torsions
$\tau R = 0.16$
,
$0.20$
and
$1.00$
, with a plane piston source of frequency
$\omega = 0.95\overline {\omega }_1^1 = 0.95 \times 1.8412/R$
, viewed from two different angles in each case. Truncation was taken at
$\alpha _{\textrm{max}} = a_{\textrm{max}} = 10$
. An animated version of this figure is available as Movie 7 in the supplementary material.
Pressure field (normalised by the Mach number) on the walls of a helical duct with torsion
$\tau R = 0.20$
, for
$M = 0.05$
(a,b), 0.10 (c,d) and 0.15 (e,f), with a plane piston source of frequency
$\omega = 0.95\overline {\omega }_1^1 = 0.95 \times 1.8412/R$
. Truncation was taken at
$\alpha _{\textrm{max}} = a_{\textrm{max}} = 10$
. An animated version of this figure is available as Movie 8 in the supplementary material.
4.8. Comparison of two and three dimensions
Width variation should not, by itself, produce significantly different effects between two and three dimensions, since the same symmetry is being broken in each case. Torsion is only possible in three dimensions, so does not invite comparison. Therefore, we will concentrate here on the difference between the effects of curvature in two and three dimensions. We consider a duct of width
$X_{\textrm{i}} = 1$
in two dimensions, and radius
$R_{\textrm{i}} = 0.5$
in three. In each case, there are straight sections of length
$7/4$
either side of a right-angled bend with bend radius
$5/4$
. We prescribe a piston source of frequency
$9/2$
: figure 14 compares the pressure fields in the linear case with those of Mach number
$0.1$
. In figure 15, we plot the spatially averaged outlet pressures against time, as a means of demonstrating the bend’s effect on a piston source. Truncation was taken at 10 spatial modes and 10 temporal modes for each calculation.
We note that while the two cases seem to steepen equally, the effective acoustic length of the bend differs between them. The three-dimensional outlet pressure is consistently ‘ahead’ of the two-dimensional outlet pressure for both linear and nonlinear regimes, suggesting that the effective bend length is longer in two dimensions. It may also be observed that the two-dimensional calculations exhibit more transverse oscillations than in three dimensions: this is unsurprising, since that is the only possible plane of oscillation in two dimensions, whereas in three dimensions there is an entire circular surface of oscillation at each
$s$
.
Forward-going pressure (normalised to the source amplitude) in an extended bend with a plane piston source of
$9/2$
, for (a) linear, two dimensions, (b) nonlinear, two dimensions, (c) linear, three dimensions and (d) nonlinear, three dimensions. Truncation was taken at
$\alpha _{\textrm{max}} = a_{\textrm{max}} = 10$
. An animated version of this figure is available as Movie 9 in the supplementary material.
Forward-going pressure (normalised to the source amplitude), averaged over the inlet, the outlet in two dimensions, and the outlet in three dimensions, in the linear regime (a) and for
$M = 0.10$
(b). Truncation was taken at
$\alpha _{\textrm{max}} = a_{\textrm{max}} = 10$
.
4.9. Effective acoustic length of a bend
We can also examine the competing effects of curvature and nonlinearity. We calculate the cross-sectional-average pressure two wavelengths along the centreline from the inlet, for a range of different bend angles and Mach numbers, while keeping the arclength along the centre of the bend constant. The boundary condition at the inlet ensures that the average pressure across the inlet will pass upwardly through zero at
$t = \pi/\omega$
, and downwardly at
$t = \pi/2\omega$
. In a straight duct, the average pressure across a surface two wavelengths along would match this, but curvature causes a deviation of these upstroke and downstroke crossings, which we call
$\Delta t$
. We then define a bend correction factor
$B := 2c\Delta t/\Delta x$
, where
$\Delta x$
is the difference in arclength between the outer wall of the duct and the centreline for that particular bend angle. This parameter is equal to 1 for a wave travelling along the outside of the bend, and −1 for the inside. We can also choose whether to use the deviation time
$\Delta t$
for the upstroke of the wave, or the downstroke. Figure 16 shows upstroke and downstroke values of
$B$
for three duct radii. Each data point on these contour plots is the result of a calculation with 10 spatial modes and 10 temporal modes.
Upstroke (a,c,e) and downstroke (b,d,f) values of
$B$
, for
$\omega X = 3$
(a,b),
$3/2$
(c,d) and
$3/4$
(e,f). Truncation was taken at
$\alpha _{\textrm{max}} = a_{\textrm{max}} = 10$
.
From the upstroke plots we see that for greater bend angles (and correspondingly tighter bends), nonlinearity plays a bigger role in forcing the wave around the outside of the bend, whereas for smaller bend angles, the bend angle takes over as the determining parameter. The downstroke plots are less conclusive, but they do exhibit the emergence of a local minimum in
$B$
as the radius is varied.
5. Conclusions
We have presented a numerically tractable mathematical framework for solving for weakly nonlinear acoustics in curved ducts of varying width, along with a number of numerical examples of its use. The framework unifies both two- and three-dimensional governing equations, and allows for torsion in three dimensions. The framework is numerically tractable in that it may be written in terms of large matrices and tensors that are invariant along the duct, unlike previous work in two dimensions where the large matrices and tensors varied with the duct geometry and therefore needed computing at each point along the duct (McTavish & Brambley Reference McTavish and Brambley2019). The method follows the multi-modal method of Félix & Pagneux (Reference Félix and Pagneux2001, Reference Félix and Pagneux2002), solving first for the admittance (a generalised ratio of acoustic velocity to acoustic pressure) throughout the duct from the outlet to the inlet, and then using this admittance to calculate the acoustic pressure and velocity from the inlet to the outlet. This has the advantage that the effects of the duct geometry are encoded in the admittance and may be analysed independently of the particular acoustic source used. The unification of two and three dimensions allows for comparison between the two cases (such as was shown in figure 14).
The method has been validated for a range of duct geometries where there are known analytical solutions or published numerical results. Nonlinear steepening is validated against the Fubini (Reference Fubini1935), Blackstock (Reference Blackstock1966) and Fay (Reference Fay1931) solutions for a straight duct in § 4.1. Width variation in the linear case is validated against the Webster (Reference Webster1919) solution for an exponential horn in § 4.4; this test case highlights the inaccuracy of approximating the full sound field by a plane wave in Webster’s solution, which can be reproduced in our framework by restricting our solution to a single spatial mode or relaxed by allowing for many more spatial modes. Curvature is validated in two dimensions in § 4.2, in the linear regime against a result of Félix & Pagneux (2001) and in the nonlinear regime against a result for the same geometry of McTavish & Brambley (Reference McTavish and Brambley2019). The three dimensional equivalent in § 4.3 is validated against a result of Félix & Pagneux (Reference Félix and Pagneux2002) in the linear regime, and against an unpublished result of McTavish (Reference McTavish2018) in the nonlinear regime. This extensive validation gives confidence that the combination of width variation, curvature, and nonlinearity is correctly modelled both in the mathematical framework and in its numerical implementation.
A number of duct geometries have then been investigated numerically. These include a demonstration of wave leakage, where the inclusion of either slight curvature or slight nonlinearity can cause waves to propagate through a duct in which they would otherwise have been totally reflected (§ 4.6). Examples with torsion in three dimensions (§ 4.7) show plane waves becoming localised on the outside of the duct curve. The unification of two- and three-dimensional equations allows for the comparison of the two in § 4.8, showing that between two and three dimensions, wave steepening appears comparable but the effective acoustic length of a bend is different.
The effective length of a bend is investigated further in two dimensions in § 4.9. By considering the time lag in a wave propagating around a bend, a bend correction factor
$B$
is introduced, taking the value
$+1$
for a wave propagating around the outside of the bend and
$-1$
for a wave propagating around the inside of the bend. These results are presented in figure 16, and suggest that for higher curvature bends nonlinearity plays a bigger role in forcing wave propagation around the outside of the bend, whereas for lower curvature bends, the curvature takes over as the determining parameter. This tool could potentially be used in the future to investigate whether resonant frequencies are likely to change when the sound amplitude is varied.
An obvious motivation of this work was the potential application to sound in brass instruments, where the brassy sound is caused by nonlinear wave steepening within the instrument. The use of the admittance (both the usual linear admittance and its weakly nonlinear extension) makes the present framework well suited to this, as a duct’s resonances could be investigated by solving for the admittance without needing to model or specify a sound source at the inlet. Moreover, the bend correction factor of § 4.9 could be used to investigate whether instruments are likely to become sharper or flatter when played at louder or quieter volumes, and potentially even whether an instrument could be designed to have a stable pitch independently of the volume at which it is played.
The numerical solution for the admittance, and then for the pressure and velocity, is relatively standard, and future work could consider other methods of numerical solution. As described in § 3.5, the infinite series of ODEs is truncated to a finite system by specifying a maximum number of temporal and spatial modes. Since this truncation prevents the cascade of energy to higher modes and leads to artificial energy accumulating in the highest non-truncated mode, an artificial viscosity is added to help dissipate energy at higher modes. This artificial viscosity is akin to the real molecular viscosity of the gas. However, the dominant energy loss mechanism for acoustics in a duct is friction with the walls, which has been previously modelled using a fractional derivative approach (see, e.g. Rendón et al. Reference Rendón, Orduña Bustamante, Narezo, Pérez-López and Sorrentini2010, and references therein); a similar technique could be investigated for inclusion in the present model. The now finite set of ODEs is then integrated here using a standard Runge–Kutta method, either with a fixed (Runge–Kutta fourth-order, RK4) or adaptive (Runge–Kutta fourth-order five-stage, RK45) method. More sophisticated numerical methods that take note of the Riccati-style nature of the equations being integrated are also possible. For example, Pagneux (Reference Pagneux2010) used a Magnus–Möbius scheme in the linear case to avoid singularities of the admittance caused by the presence of pressure nodes in the duct. Since the extension of Magnus schemes to nonlinear equations is possible (Casas & Iserles Reference Casas and Iserles2006), potential future research might investigate the use of a weakly nonlinear Magnus–Möbius scheme for at least the admittance integration.
In order to calculate the admittance within the duct, a known admittance must be applied, for example at the duct outlet. In all of the examples presented here, a non-reflecting admittance representing an infinite straight duct of constant width was used as the outlet admittance, as discussed in § 3.2. One obvious outlet admittance for a musical instrument would be that of a duct exiting into an infinite space; such a situation may be modelled in the linear case either exactly using a Wiener–Hopf technique (Munt Reference Munt1977) or approximately using an outer duct with absorbing walls (e.g. Félix et al. Reference Félix, Doc and Boucher2018), but the authors are not aware of a nonlinear, or even weakly nonlinear, equivalent. Both the open duct exit admittance and the directivity pattern of the far-field radiation would be interesting avenues of future research.
Other possibilities for future research might involve the inclusion of a mean flow within the duct (e.g. Mangin et al. Reference Mangin, Daroukh and Gabard2023), ducts with acoustically lined walls (e.g. Bi et al. Reference Bi, Pagneux, Lafarge and Aurégan2007), the investigation of other nonlinear effects such as resonant triads (e.g. Protas et al. Reference Protas, Kang and Bustamante2024) and musical dynamics owing to the feedback loop between a musician’s lips and the reflections from the waves that they generate, or other modal representations that might yield faster numerical convergence (e.g. Maurel, Mercier & Pagneux Reference Maurel, Mercier and Pagneux2014).
Supplementary material and movies
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2026.11220.
Acknowledgements
For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.
Funding
F.J. was funded through the Warwick Mathematics Institute Centre for Doctoral Training, and gratefully acknowledges the support of the University of Warwick and the UK Engineering and Physical Sciences Research Council (EPSRC grant EP/W523793/1). E.J.B. gratefully acknowledges the support of the UK Engineering and Physical Sciences Research Council (EPSRC grant EP/V002929/1).
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Spatial projection in more detail
A.1. Two-dimensional spatial projection
The goal here is to project (2.24) onto the basis of spatial modes (2.41). This is achieved by multiplying each equation by a mode
$\psi _\alpha$
and integrating over a duct cross-section. If the expansions (2.37) are then employed, the acoustic mode coefficients may then be factored out of the integral, which becomes a matrix with rows in
$\alpha$
and columns in the expansion dummy variable
$\beta$
, dependent only on the local geometry. However, the Neumann boundary condition (2.40) prevents us from expanding any first derivatives of acoustical quantities on the boundary, which means we may not expand second derivatives of acoustical quantities within integrals, and must instead use integration by parts to remove the
$x$
-derivatives and then expand once it is safe to so.
If we start with the left-hand side of (2.24a), we get
\begin{align} \int _{X_-}^{X_+}\psi _\alpha \left \{\frac {\partial U^a}{\partial s} - {\textrm{i}}\kern0.5pt a\omega \left [h_s\left (1 + \frac {1}{a^2\omega ^2}\frac {\partial ^2}{\partial x^2}\right ) - \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial x}\right ]P^a\right \}{\rm d} x = \underbrace {\frac {{\rm d} }{{\rm d} s}\left (\int _{X_-}^{X_+}\psi _\alpha U^a\,{\rm d} x\right )}_{\fbox{$1$}}\nonumber \\ - \underbrace {\left [X_+'\left (\psi _\alpha U^a\right )\bigg |_{X_+} - X_-'\left (\psi _\alpha U^a\right )\bigg |_{X_-}\right ]}_{\fbox{$2$}} - \underbrace {\frac {{\textrm{i}}\kern0.5pt a\omega }{a^2\omega ^2}\left [\psi _\alpha h_s\frac {\partial P^a}{\partial x} - \frac {\partial (\psi _\alpha h_s)}{\partial x}P^a - \kappa \psi _\alpha P^a\right ]_{X_-}^{X_+}}_{\fbox{$3$}} \nonumber \\ - \int _{X_-}^{X_+}\underbrace {\frac {\partial \psi _\alpha }{\partial s}U^a}_{\fbox{$4$}} + \underbrace {{\textrm{i}}\kern0.5pt a\omega \left (h_s\psi _\alpha + \frac {1}{a^2\omega ^2}\frac {\partial ^2(\psi _\alpha h_s)}{\partial x^2} + \frac {\kappa }{a^2\omega ^2}\frac {\partial \psi _\alpha }{\partial x}\right )P^a}_{\fbox{$5$}}{\rm d} x. \end{align}
Terms
$\fbox{$1$}$
,
$\fbox{$2$}$
and
$\fbox{$4$}$
are the result of bringing the partial
$s$
-derivative on the
$U^a$
outside of the integral, whereas terms
$\fbox{$3$}$
and
$\fbox{$5$}$
are the result of removing all of the
$x$
-derivatives from the
$P^a$
via integration by parts (in particular, the second derivative has been removed by a double application of integration-by-parts to produce the first two constituent terms of
$\fbox{$3$}$
). Term
$\fbox{$1$}$
is an exact
$s$
-derivative of the projection of
$U^a$
onto the basis, i.e.
\begin{equation} \frac {{\rm d} }{{\rm d} s}\left (\int _{X_-}^{X_+}\psi _\alpha U^a{\rm d} x\right ) = \frac {{\rm d} }{{\rm d} s}\left (\int _{X_-}^{X_+}\psi _\alpha \sum _\beta U_\beta ^a\psi _\beta {\rm d} x\right ) = \frac {{\rm d} }{{\rm d} s}\left (\sum _\beta \delta _{\alpha \beta }U_\beta ^a\right ) = \frac {{\rm d} U_\alpha ^a}{{\rm d} s}, \end{equation}
by mode orthogonality. Term
$\fbox{$2$}$
combines with the first term of
$\fbox{$3$}$
to form the left-hand side of the hard-walled boundary condition (2.25), leaving a term that is
$O(M^2)$
, while the second two terms of
$\fbox{$3$}$
cancel with one another due to the Neumann condition on the duct modes. Finally,
$\fbox{$5$}$
’s second term, by the Leibniz rule, makes -2 lots of the
$\fbox{$5$}$
third term added to a term involving the second derivative of a spatial mode, which by (2.38) may be turned into multiplication by a squared eigenvalue. Thus, expanding into spatial modes, we have
\begin{align} \int _{X_-}^{X_+}\psi _\alpha &\left \{\frac {\partial U^a}{\partial s} - {\textrm{i}}\kern0.5pt a\omega \left [h_s\left (1 + \frac {1}{a^2\omega ^2}\frac {\partial ^2}{\partial x^2}\right ) - \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial x}\right ]P^a\right \}{\rm d} x \notag \\ =& \sum _\beta \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s} - \int _{X_-}^{X_+}\frac {\partial \psi _\alpha }{\partial s}\psi _\beta {\rm d} x\right )U_\beta ^a + \frac {1}{{\textrm{i}}\kern0.5pt a\omega }\left [h_s\psi _\alpha \frac {\partial Q^a}{\partial x}\right ]_{X_-}^{X_+} \notag \\ &- {\textrm{i}}\kern0.5pt a\omega \sum _\beta \left [\left (1 - \frac {\lambda _\alpha ^2}{a^2\omega ^2X^2}\right )\int _{X_-}^{X_+}h_s\psi _\alpha \psi _\beta\, {\rm d} x - \frac {\kappa }{a^2\omega ^2}\int _{X_-}^{X_+}\frac {\partial \psi _\alpha }{\partial x}\psi _\beta\, {\rm d} x\right ]P_\beta ^a. \end{align}
Apart from the term involving
$Q^a$
(which will be dealt with later), we have integrals, with two subscripts, summing over their second subscripts with acoustical quantities that have only one subscript. In other words, this is standard matrix-vector multiplication, so we may drop the sigmas in employment of the summation convention (for clarity, in all of the expressions henceforth
$\alpha$
is the free suffix, while,
$\beta$
,
$\gamma$
and
$\delta$
are dummy suffices for use in summation). Furthermore, we will define a shorthand for any mode-integral matrices: the integral of two modes is written as
$\varPsi$
, with Greek letter subscripts denoting the mode numbers, while bracketed subscripts denote derivatives of the modes in question (curly brackets are
$s$
-derivatives and square brackets are
$x$
-derivatives, with a square bracket around all subscripts meaning that this is a boundary term). Finally, any non-mode functions present in the integral appear in square brackets after the subscripts. This notation is the same as was used in McTavish & Brambley (Reference McTavish and Brambley2019). For example,
\begin{align} &\varPsi _{\{\alpha \}[\beta ]}[f(s,x)] \quad \text{would mean}\quad \int _{X_-}^{X_+}\frac {\partial \psi _\alpha }{\partial s}\frac {\partial \psi _\beta }{\partial x}f(s,x){\rm d} x. \end{align}
With this notation in hand, the projected left-hand side becomes
\begin{align} \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s} - \varPsi _{\{\alpha \}\beta }\right )U_\beta ^a & + \frac {1}{{\textrm{i}}\kern0.5pt a\omega }\left [h_s\psi _\alpha \frac {\partial Q^a}{\partial x}\right ]_{X_-}^{X_+} \nonumber \\& - {\textrm{i}}\kern0.5pt a\omega \left [\left (\delta _{\alpha \gamma } - \frac {\boldsymbol{\varLambda }_{\alpha \gamma }^2}{a^2\omega ^2X^2}\right )\varPsi _{\gamma \beta }[h_s] - \frac {\kappa }{a^2\omega ^2}\varPsi _{[\alpha ]\beta }\right ]P_\beta ^b = O(M^2), \end{align}
where
$\boldsymbol{\varLambda }$
is a matrix with entries given by
$\boldsymbol{\varLambda }_{\alpha \gamma } = \lambda _{\alpha }\delta _{\alpha \gamma }$
(no sum).
We now turn to the right-hand side of (2.24a ). Its first term projects very easily, giving
\begin{equation} \int _{X_-}^{X_+}-{\textrm{i}}\kern0.5pt a\omega \beta_0 h_s\psi _\alpha \sum _bP^{a-b}P^b{\rm d} x = -{\textrm{i}}\kern0.5pt a\omega \beta_0\varPsi _{\alpha \beta \gamma }[h_s]\sum _bP_\beta ^{a - b}P_\gamma ^b, \end{equation}
where the
$\varPsi$
matrix notation has been naturally extended to third-rank tensors, which contract with acoustical quantities on subscripts 2 and 3. The next term on the right-hand side becomes
\begin{align}& {\textrm{i}}\kern0.5pt a\omega \int _{X_-}^{X_+}\psi _\alpha \left [h_s\left (1 - \frac {1}{a^2\omega ^2}\frac {\partial ^2}{\partial x^2}\right ) + \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial x}\right ]Q^a{\rm d} x \nonumber \\[4pt]& \qquad \qquad = \frac {{\textrm{i}}\kern0.5pt a\omega }{a^2\omega ^2}\bigg [-\psi _\alpha h_s\frac {\partial Q^a}{\partial x} + \frac {\partial (\psi _\alpha h_s)}{\partial x}Q^a + \kappa \psi _\alpha Q^a\bigg ]_{X_-}^{X_+} \nonumber \\[4pt]& \qquad \qquad \quad + {\textrm{i}}\kern0.5pt a\omega \int _{X_-}^{X_+}Q^a\left [\left (1 - \frac {1}{a^2\omega ^2}\frac {\partial ^2}{\partial x^2}\right )(h_s\psi _\alpha ) - \frac {\kappa }{a^2\omega ^2}\frac {\partial \psi _\alpha }{\partial x}\right ]{\rm d} x. \end{align}
This is more or less the same manipulation as took place on the left-hand side, only now with
$Q^a$
in the place of
$P^a$
and some sign changes. After the same cancellations have been made, and we switch to the
$\varPsi$
notation (which requires the terms making up
$Q^a$
to be expanded into spatial modes), we get
\begin{align} {\textrm{i}}\kern0.5pt a\omega \int _{X_-}^{X_+}\psi _\alpha \left [h_s\left (1 - \frac {1}{a^2\omega ^2}\frac {\partial ^2}{\partial x^2}\right ) + \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial x}\right ]Q^a\,{\rm d} x = \frac {1}{{\textrm{i}}\kern0.5pt a\omega }\left [h_s\psi _\alpha \frac {\partial Q^a}{\partial x}\right ]_{X_-}^{X_+} \nonumber \\ + {\textrm{i}}\kern0.5pt a\omega \bigg [\left (\delta _{\alpha \delta } + \frac {\boldsymbol{\varLambda }_{\alpha \delta }^2}{a^2\omega ^2X^2}\right )\varPsi _{\delta \beta \gamma }[h_s] + \frac {\kappa }{a^2\omega ^2}\varPsi _{[\alpha ]\beta \gamma }\bigg ]\frac {1}{2}\sum _bP_\beta ^{a - b}P_\gamma ^b - U_\beta ^{a - b}U_\gamma ^b \nonumber \\ + {\textrm{i}}\kern0.5pt a\omega \bigg [\left (\delta _{\alpha \delta } + \frac {\boldsymbol{\varLambda }_{\alpha \delta }^2}{a^2\omega ^2X^2}\right )\varPsi _{\delta [\beta ][\gamma ]}[h_s] + \frac {\kappa }{a^2\omega ^2}\varPsi _{[\alpha ][\beta ][\gamma ]}\bigg ]\sum _b\frac {P_\beta ^{a - b}P_\gamma ^b}{2(a - b)b\omega ^2}. \end{align}
The first term here is exactly the leftover boundary term from the left-hand side; this will therefore cancel when they are combined, meaning the projected mass conservation equation is
\begin{align} \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s}\right .&\left .{}\!\!- \varPsi _{\{\alpha \}\beta }\right )U_\beta ^a - {\textrm{i}}\kern0.5pt a\omega \left [\left (1 - \frac {\lambda _\alpha ^2}{a^2\omega ^2X^2}\right )\varPsi _{\alpha \beta }[h_s] - \frac {\kappa }{a^2\omega ^2}\varPsi _{[\alpha ]\beta }\right ]P_\beta ^b \nonumber\\ =&\, {\textrm{i}}\kern0.5pt a\omega \sum _b\Bigg \{- \beta_0\varPsi _{\alpha \beta \gamma }[h_s]P_\beta ^{a - b}P_\gamma ^b \nonumber\\ &+ \bigg [\left (\delta _{\alpha \delta } + \frac {\boldsymbol{\varLambda }_{\alpha \delta }^2}{a^2\omega ^2X^2}\right )\varPsi _{\delta \beta \gamma }[h_s] + \frac {\kappa }{a^2\omega ^2}\varPsi _{[\alpha ]\beta \gamma }\bigg ]\frac {P_\beta ^{a - b}P_\gamma ^b - U_\beta ^{a - b}U_\gamma ^b}{2} \nonumber\\ &+ \bigg [\left (\delta _{\alpha \delta } + \frac {\boldsymbol{\varLambda }_{\alpha \delta }^2}{a^2\omega ^2X^2}\right )\varPsi _{\delta [\beta ][\gamma ]}[h_s] + \frac {\kappa }{a^2\omega ^2}\varPsi _{[\alpha ][\beta ][\gamma ]}\bigg ]\frac {P_\beta ^{a - b}P_\gamma ^b}{2(a - b)b\omega ^2}\Bigg \}. \end{align}
Projecting (2.24b) proves easier, since there are no boundary-term cancellations; care is only needed when dealing with the second
$x$
-derivative on the right-hand side. The left-hand side requires no integration-by-parts this time, due to the absence of any
$x$
-derivatives, so we simply have
\begin{equation} \int _{X_-}^{X_+}\psi _\alpha \left (\frac {\partial P^a}{\partial s} - {\textrm{i}}\kern0.5pt a\omega h_sU^a\right ){\rm d} x = \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s} + \varPsi _{\alpha \{\beta \}}\right )P_\beta ^a - {\textrm{i}}\kern0.5pt a\omega \varPsi _{\alpha \beta }[h_s]U_\beta ^a. \end{equation}
The term involving
$(a - b)$
on the right-hand side is dealt with very easily, and the final term is effectively a repeat of the final term in (A8). This leaves only the term with the
$\partial ^2/\partial x^2$
, which must be integrated by parts before an expansion of the pressure into spatial modes may take place
\begin{align} \int _{X_-}^{X_+}{\textrm{i}}\kern0.5pt\omega h_s\psi _\alpha \sum _b\left [-bU^{a - b}\frac {1}{b^2\omega ^2}\frac {\partial ^2P^b}{\partial x^2}\right ]{\rm d} x = -{\textrm{i}}\kern0.5pt\omega \left [h_s\psi _\alpha \sum _b\frac {b}{b^2\omega ^2}U^{a - b}\frac {\partial P^b}{\partial x}\right ]_{X_-}^{X_+} \nonumber \\ + \int _{X_-}^{X_+} {\textrm{i}}\kern0.5pt\omega \sum _b\left [\frac {b}{b^2\omega ^2}\frac {\partial \left (h_s\psi _\alpha U^{a - b})\right )}{\partial x}\sum _\gamma P_\gamma ^b\frac {\partial \psi _\gamma }{\partial x}\right ]{\rm d} x. \end{align}
We know from the hard-walled boundary condition that the
$\partial P^b/\partial x$
in the boundary term must turn into a
$U^b$
(i.e. something expansible), and having now (legally) expanded the
$P^b$
inside the integral, we may now simply reverse the integration-by-parts (and employ the summation convention over
$\gamma$
), giving
\begin{align} \int _{X_-}^{X_+}{\textrm{i}}\kern0.5pt\omega h_s\psi _\alpha \sum _b\left (-bU^{a - b}\frac {1}{b^2\omega ^2}\frac {\partial ^2P^b}{\partial x^2}\right ){\rm d} x = \left [X_+\varPsi _{[\alpha \beta \gamma ]}^+ - X_-\varPsi _{[\alpha \beta \gamma ]}^-\right ]\sum _bU_\beta ^{a - b}U_\gamma ^b \nonumber \\ + {\textrm{i}}\kern0.5pt\omega \sum _b\left [\frac {b}{b^2\omega ^2}h_s\psi _\alpha U^{a - b}P_\gamma ^b\frac {\partial \psi _\gamma }{\partial x}\right ]_{X_-}^{X_+} + {\textrm{i}}\kern0.5pt\omega \sum _b\frac {b\boldsymbol{\varLambda }_{\delta \gamma }^2}{b^2\omega ^2X^2}\varPsi _{\alpha \beta \delta }[h_s]U_\beta ^{a - b}P_\gamma ^b. \end{align}
The
$\partial \psi _\gamma /\partial x$
in the second boundary term causes it to vanish, by the Neumann condition on spatial modes. Thus, this term has been dealt with, and the projected momentum equation is
\begin{align}& \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s} + \varPsi _{\alpha \{\beta \}}\right )P_\beta ^a - {\textrm{i}}\kern0.5pt a\omega \varPsi _{\alpha \beta }[h_s]U_\beta ^a = \left [X_+\varPsi _{[\alpha \beta \gamma ]}^+ - X_-\varPsi _{[\alpha \beta \gamma ]}^-\right ]\sum _bU_\beta ^{a - b}U_\gamma ^b \nonumber \\[3pt]& + {\textrm{i}}\kern0.5pt\omega \sum _b\left \{\varPsi _{\alpha \beta \delta }[h_s]\left [(a - b)\delta _{\delta \gamma } - b\left (\delta _{\delta \gamma } - \frac {\boldsymbol{\varLambda }_{\delta \gamma }^2}{b^2\omega ^2X^2}\right )\right ] + \frac {1}{b\omega ^2}\varPsi _{\alpha [\beta ][\gamma ]}[h_s]\right \} \!U_\beta ^{a - b}P_\gamma ^b, \end{align}
where the ‘
$+$
’ and ‘
$-$
’ superscripts on boundary terms indicate that only a single side of the boundary is being evaluated.
In theory, we have now done all of the work required to turn these equations from PDEs to ODEs and the next step would be to put them into a computer and solve them. However, more work may now be done to these equations in order to make them numerically efficient to solve. At present, we have a countably infinite set of coupled vector ODEs in
$s$
with matricial coefficients that also depend on
$s$
, meaning that a solver would need to update the value of each matrix at each point along the duct. This is missing a trick, however, since each of these matrices’
$s$
-dependence may be factored out into the form of
$s$
-dependent scalars such as
$\kappa$
, while the matrices themselves are constant, and thus only need to be defined before solving. This also has the advantage of added clarity, since it will clearly highlight which terms in the equations are responsible for different geometrical irregularities. In two dimensions, we have
$s$
-dependence in the integration limits
$X_-$
and
$X_+$
, which is removed via the substitution
$\xi = (x - X_-)/X$
, so that the integration limits are now 0 and 1. Expanding
$h_s$
and factoring the
$\kappa$
out then just leaves
$1/\sqrt {X}$
in front of the modes, which can also be factored out. If we apply this technique to one of the matrices from (A13), we get
\begin{align} \varPsi _{\alpha \beta }[h_s] &= \int _{X_-}^{X_+}\frac {C_\alpha C_\beta }{X}h_s\cos \left (\lambda _\alpha \frac {x - X_-}{X}\right )\cos \left (\lambda _\beta \frac {x - X_-}{X}\right ){\rm d} x \nonumber \\[3pt]& = (1 - \kappa X_-)\delta _{\alpha \beta } - \kappa X\int _0^1C_\alpha C_\beta \xi \cos (\lambda _\alpha \xi )\cos (\lambda _\beta \xi ){\rm d} \xi \nonumber \\[3pt] &= (1 - \kappa X_-)\delta _{\alpha \beta } - \kappa X\varXi _{\alpha \beta }[\xi ] ,\end{align}
where we have adopted a new notation
$\varXi$
for the constant matrices, analogous to the
$\varPsi$
notation but now with different limits and a different integration variable. In order to deal with
$s$
-derivatives, we note that
\begin{align} \frac {\partial \psi _\alpha }{\partial s} = -\frac {X'}{2X}\psi _\alpha - \frac {X'}{X}(x - X_-)\frac {\partial \psi _\alpha }{\partial x} - X_-'\frac {\partial \psi _\alpha }{\partial x} = - \frac {X'}{2X}\left (\psi _\alpha + 2\xi \frac {\partial \psi _\alpha }{\partial \xi }\right ) - \frac {X_-'}{X}\frac {\partial \psi _\alpha }{\partial \xi },\\[-18pt]\nonumber \end{align}
and a more convenient form may be found for some of the right-hand-side terms involving derivatives on multiple indices, by use of the Leibniz rule. The two relations necessary for this are
\begin{align} \varPsi _{\alpha [\beta ][\gamma ]}[h_s] &= \frac {\lambda _\beta ^2 + \lambda _\gamma ^2 - \lambda _\alpha ^2}{2X^2}\varPsi _{\alpha \beta \gamma }[h_s] - \kappa \varPsi _{[\alpha ]\beta \gamma } + \frac {\kappa }{2}\varPsi _{[\alpha \beta \gamma ]}, \\[-12pt] \nonumber \end{align}
\begin{align} \varPsi _{[\alpha ][\beta ][\gamma ]} &= \frac {\lambda _{\beta }^2 + \lambda _\gamma ^2 - \lambda _\alpha ^2}{2X^2}\varPsi _{[\alpha ]\beta \gamma } + \frac {\lambda _\alpha ^2}{2X^2}\varPsi _{[\alpha \beta \gamma ]}. \\[9pt] \nonumber \end{align}
In § 2.6 we introduced a compact notation which computationally and algebraically simplified the equations. In the new notation, (A9) and (A13) can be rewritten as
\begin{align} & \Bigg [\frac {{\rm d} }{{\rm d} s} + \frac {X'}{2X} \unicode{x1D652} + \frac {X_-'}{X}\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}\Bigg ]\boldsymbol{u}^a - {\textrm{i}}\kern0.5pt a\omega \Bigg [\!\!\left (\! \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X^2}\!\right )\!\big((1 - \kappa X_-) \unicode{x1D644} - \kappa X \unicode{x1D63C}\big) - \frac {\kappa \kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}}{a^2\omega ^2X} \!\Bigg ]\boldsymbol{p}^a \notag \\ =& \frac {{\textrm{i}}\kern0.5pt a\omega }{\sqrt {X}}\sum _b\Bigg \{- \beta_0((1 - \kappa X_-)\mathcal{I} - \kappa X\mathcal{A})\langle \boldsymbol{p}^{a - b},\boldsymbol{p}^b\rangle \notag \\& \quad+\Bigg [\left ( \unicode{x1D644} + \frac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X^2}\right ) \big((1 - \kappa X_-)\mathcal{I} - \kappa X\mathcal{A}\big) + \frac {\kappa \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}}{a^2\omega ^2X}\Bigg ]\frac {\langle \boldsymbol{p}^{a - b},\boldsymbol{p}^b\rangle - \langle \boldsymbol{u}^{a - b},\boldsymbol{u}^b\rangle }{2} \notag \\& \quad + \left [\frac {\left ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X^2}\right ) \big((1 - \kappa X_-)\mathcal{I}^\lambda - \kappa X\mathcal{A}^\lambda\big) + \dfrac {\kappa \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda }{a^2\omega ^2X}}{2(a - b)b\omega ^2X^2} \right ]\langle \boldsymbol{p}^{a - b},\boldsymbol{p}^b\rangle \Bigg \}, \\[-12pt] \nonumber \end{align}
\begin{align} & \Bigg [\frac {{\rm d} }{{\rm d} s} - \frac {X'}{2X} \unicode{x1D652}^{{\,T}} - \frac {X_-'}{X}\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}^{{T}}\Bigg ]\boldsymbol{p}^a - {\textrm{i}}\kern0.5pt a\omega \big [(1 - \kappa X_-) \unicode{x1D644} - \kappa X \unicode{x1D63C}\big ]\boldsymbol{u}^a \notag \\\quad=& \frac {1}{\sqrt {X}}\sum _b\Bigg \{\bigg (\frac {X_+'}{X}\overline {\mathcal{W}}^+ - \frac {X_-'}{X}\overline {\mathcal{W}}^-\bigg )\langle \boldsymbol{u}^{a - b},\boldsymbol{u}^b\rangle \notag \\& \quad+ {\textrm{i}}\kern0.5pt\omega \Bigg [\big ((1 - \kappa X_-)\mathcal{I} - \kappa X\mathcal{A}\big ) \left \langle \unicode{x1D644},(a - b) \unicode{x1D644} - b\left ( \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{b^2\omega ^2X^2}\right )\right \rangle \notag \\&\quad\qquad \quad + \frac {(1 - \kappa X_-)\mathcal{I}^\lambda - \kappa X\mathcal{A}^\lambda }{b\omega ^2X^2}\Bigg ]\langle \boldsymbol{u}^{a - b},\boldsymbol{p}^b\rangle \Bigg \}, \\[9pt] \nonumber \end{align}
where
\begin{align} \unicode{x1D644} &= \delta _{\alpha \beta }, & \unicode{x1D652}_{\alpha \beta } &= \delta _{\alpha \beta } + 2\varXi _{[\alpha ]\beta }[\xi ], \\[-12pt] \nonumber \end{align}
\begin{align} \unicode{x1D63C}_{\alpha \beta } &= \varXi _{\alpha \beta }[\xi ], & \kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}_{\alpha \beta } &= \varXi _{[\alpha ]\beta }, \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{I}_{\alpha \beta \gamma } &= \varXi _{\alpha \beta \gamma }, & \mathcal{A}_{\alpha \beta \gamma } &= \varXi _{\alpha \beta \gamma }[\xi ], \\[-12pt] \nonumber \end{align}
\begin{align} \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}_{\alpha \beta \gamma } &= \varXi _{[\alpha ]\beta \gamma }, & \kern2.5pt \overline {\kern-2.5pt \mathcal{A}}_{\alpha \beta \gamma } &= \varXi _{[\alpha \beta \gamma ]}, \\[-12pt] \nonumber \end{align}
\begin{align} \overline {\mathcal{W}}_{\alpha \beta \gamma }^+ &= \varXi _{[\alpha \beta \gamma ]}^+, & \overline {\mathcal{W}}_{\alpha \beta \gamma }^- &= \varXi _{[\alpha \beta \gamma ]}^-, \\[-6pt] \nonumber \end{align}
\begin{align} \mathcal{I}^\lambda &= \mathcal{I}\frac {\langle \boldsymbol{\varLambda }^2, \unicode{x1D644}\,\rangle + \langle \unicode{x1D644},\boldsymbol{\varLambda }^2\rangle - \boldsymbol{\varLambda }^2}{2}, \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{A}^\lambda &= \mathcal{A}\frac {\langle \boldsymbol{\varLambda }^2, \unicode{x1D644}\,\rangle + \langle \unicode{x1D644},\boldsymbol{\varLambda }^2\rangle - \boldsymbol{\varLambda }^2}{2} + \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}} - \frac {\kern2.5pt \overline {\kern-2.5pt \mathcal{A}}}{2}, \\[-12pt] \nonumber \end{align}
\begin{align} \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda &= \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}\frac {\langle \boldsymbol{\varLambda }^2, \unicode{x1D644}\,\rangle + \langle \unicode{x1D644},\boldsymbol{\varLambda }^2\rangle - \boldsymbol{\varLambda }^2}{2} + \frac {\boldsymbol{\varLambda }^2}{2}\kern2.5pt \overline {\kern-2.5pt \mathcal{A}}. \\[9pt] \nonumber \end{align}
The names of these matrices have been chosen to represent the geometric irregularities whose influence upon the equations they modulate, i.e.
$ \unicode{x1D652}$
for ‘width’, and
$ \unicode{x1D63C}$
for ‘annularity’ (one matrix,
$\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}$
is used for both, so a single name had to be picked). Calligraphic tensors are related to sans-serif matrices with the same name, i.e.
$\mathcal{I}_{\alpha \beta 0} = \unicode{x1D644}_{\alpha \beta }$
. (A17) are (2.48) in the main text.
A.2. Three-dimensional spatial projection
Projection of the governing equations is complicated by the extra dimension, however, fortunately, all boundary terms in
$\theta$
disappear due to periodicity of the spatial modes. Thus, we only need to worry about boundary terms in
$r$
. To help with this, some quick formulae are introduced for integration by parts (over a duct cross-section
$S$
, with differential
${\rm d} S = r{\rm d} r{\rm d} \theta$
) involving the transverse gradient and FSND
\begin{gather} \begin{aligned} \iint _Sf{\nabla }_{{t}}^2g\,{\rm d} S &= \int _0^{2 \pi}\left [rf\frac {\partial g}{\partial r}\right ]_0^R{\rm d} \theta - \iint _S\boldsymbol{\nabla }_{{t}}f\boldsymbol{\cdot }\boldsymbol{\nabla }_{{t}}g\,{\rm d} S \\ &= \int _0^{2 \pi}\left [r\left (f\frac {\partial g}{\partial r} - \frac {\partial f}{\partial r}g\right )\right ]_0^R{\rm d} \theta + \iint _Sg{\nabla }_{{t}}^2f\,{\rm d} S, \end{aligned} \\[-28pt] \nonumber \end{gather}
\begin{gather} \iint _Sf\frac {\partial g}{\partial \boldsymbol{n}}{\rm d} S = \int _0^{2 \pi}\bigg [r\cos \phi fg\bigg ]_0^R{\rm d} \theta - \iint _Sg\frac {\partial f}{\partial \boldsymbol{n}}\,{\rm d} S. \\[-2pt] \nonumber \end{gather}
Note that we introduce here a shorthand
$\phi$
for use in this appendix, defined as
$\phi := \theta - \theta _0(s)$
. Using the formulae above, we can project the left-hand side of (2.35a)
\begin{align} \iint _S\psi _\alpha \left \{\frac {\partial U^a}{\partial s} - {\textrm{i}}\kern0.5pt a\omega \left [h_s\left (1 + \frac {{\nabla }_{{t}}^2}{a^2\omega ^2}\right ) - \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial \boldsymbol{n}}\right ]P^a\right \}{\rm d} S = \underbrace {\frac {{\rm d} }{{\rm d} s}\left (\iint \psi _\alpha U^a{\rm d} S\right )}_{\fbox{$1$}} \nonumber \\ - \underbrace {\int _0^{2 \pi}R'\left (r\psi _\alpha U^a\right )\bigg |_R}_{\fbox{$2$}} + \underbrace {\frac {{\textrm{i}}\kern0.5pt a\omega }{a^2\omega ^2}\left [r\psi _\alpha h_s\frac {\partial P^a}{\partial r} - r\frac {\partial (\psi _\alpha h_s)}{\partial r}P^a - \kappa r\cos \phi \psi _\alpha P^a\right ]_0^R}_{\fbox{$3$}}{\rm d} \theta \nonumber \\ -\iint _S\underbrace {\frac {\partial \psi _\alpha }{\partial s}U^a}_{\fbox{$4$}} + \underbrace {{\textrm{i}}\kern0.5pt a\omega \left (h_s\psi _\alpha + \frac {{\nabla }_{{t}}^2}{a^2\omega ^2}(h_s\psi _\alpha ) + \frac {\kappa }{a^2\omega ^2}\frac {\partial \psi _\alpha }{\partial \boldsymbol{n}}\right )P^a}_{\fbox{$5$}}{\rm d} S. \end{align}
As in (A1), we have brought the
$s$
-derivative on the
$U^a$
outside of the integral, and removed all of the transverse derivatives from the pressure via integration by parts, so that all derivatives now act on the
$\psi _\alpha$
mode. Once again term
$\fbox{$2$}$
combines with the first bit of term
$\fbox{$3$}$
to form the linear part of the boundary condition (2.36), while the rest of term
$\fbox{$3$}$
cancels, and term
$\fbox{$5$}$
changes after an application of the Leibniz rule, leaving
\begin{align} \iint _S\psi _\alpha &\left \{\frac {\partial U^a}{\partial s} - {\textrm{i}}\kern0.5pt a\omega \left [h_s\left (1 + \frac {{\nabla }_{{t}}^2}{a^2\omega ^2}\right ) - \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial \boldsymbol{n}}\right ]P^a\right \}{\rm d} S \notag \\ =& \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s} - \iint _S\frac {\partial \psi _\alpha }{\partial s}\psi _\beta {\rm d} S\right )U_\beta ^a + \frac {1}{{\textrm{i}}\kern0.5pt a\omega }\left [rh_s\psi _\alpha \frac {\partial Q^a}{\partial r}\right ]_0^R \notag \\ &- {\textrm{i}}\kern0.5pt a\omega \left [\left (1 - \frac {\lambda _\alpha ^2}{a^2\omega ^2R^2}\right )\iint _Sh_s\psi _\alpha \psi _\beta\, {\rm d} S - \frac {\kappa }{a^2\omega ^2}\iint _S\frac {\partial \psi _\alpha }{\partial \boldsymbol{n}}\psi _\beta\, {\rm d} S\right ]P_\beta ^a. \end{align}
We can use the same shorthand for these matrix integrals as we did earlier, introducing ‘
$|\bullet |$
’ as a subscript notation for the FSND, and round brackets for
$\theta$
-derivatives, so that the projected left-hand side of (2.35a
) becomes
\begin{align} & \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s} - \varPsi _{\{\alpha \}\beta }[r]\right )U_\beta ^a + \frac {1}{{\textrm{i}}\kern0.5pt a\omega }\int _0^{2 \pi}\left [rh_s\psi _\alpha \frac {\partial Q^a}{\partial r}\right ]_0^R{\rm d} \theta \nonumber \\& \qquad \qquad - {\textrm{i}}\kern0.5pt a\omega \left [\left (\delta _{\alpha \gamma } - \frac {\boldsymbol{\varLambda }_{\alpha \gamma }^2}{a^2\omega ^2R^2}\right )\varPsi _{\gamma \beta }[rh_s] - \frac {\kappa }{a^2\omega ^2}\varPsi _{|\alpha |\beta }[r]\right ]P_\beta ^a = O(M^2) .\end{align}
Turning now to the right-hand side of (2.35a
), we have an easy first term, followed by some derivatives of
$Q^a$
, which will be removed by integration-by-parts as they just were with the pressure on the left-hand side
\begin{align} {\textrm{i}}\kern0.5pt a\omega \iint _S\psi _\alpha & \left [h_s\left (1 - \frac {{\nabla }_{{t}}^2}{a^2\omega ^2}\right ) + \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial \boldsymbol{n}}\right ]Q^a{\rm d} S \nonumber \\& = \int _0^{2 \pi}\frac {{\textrm{i}}\kern0.5pt a\omega }{a^2\omega ^2}\bigg [-rh_s\psi _\alpha \frac {\partial Q^a}{\partial r} + r\frac {\partial (h_s\psi _\alpha )}{\partial r}Q^a + \kappa r\cos \phi \psi _\alpha Q^a\bigg ]_0^R{\rm d} \theta \nonumber \\& \qquad \qquad+ {\textrm{i}}\kern0.5pt a\omega \iint _S Q^a\left [\left (1 - \frac {{\nabla }_{{t}}^2}{a^2\omega ^2}\right )(h_s\psi _\alpha ) - \frac {\kappa }{a^2\omega ^2}\frac {\partial \psi _\alpha }{\partial \boldsymbol{n}}\right ]{\rm d} S .\end{align}
Switching to the
$\varPsi$
notation and making some cancellations, we get
\begin{align} & {\textrm{i}}\kern0.5pt a\omega \iint _S\psi _\alpha \bigg [h_s\left (1 - \frac {{\nabla }_{{t}}^2}{a^2\omega ^2}\right ) + \frac {\kappa }{a^2\omega ^2}\frac {\partial }{\partial \boldsymbol{n}}\bigg ]Q^a{\rm d} S = \frac {1}{{\textrm{i}}\kern0.5pt a\omega }\int _0^{2 \pi}\left [rh_s\psi _\alpha \frac {\partial Q^a}{\partial r}\right ]_0^R{\rm d} \theta \nonumber \\& \quad + {\textrm{i}}\kern0.5pt a\omega \Bigg \{\bigg [\left (\delta _{\alpha \delta } + \frac {\boldsymbol{\varLambda }_{\alpha \delta }^2}{a^2\omega ^2R^2}\right )\varPsi _{\delta \beta \gamma }[rh_s] + \frac {\kappa }{a^2\omega ^2}\varPsi _{|\alpha |\beta \gamma }[r]\bigg ]\sum _b\frac {P_\beta ^{a - b}P_\gamma ^b - U_\beta ^{a - b}U_\gamma ^b}{2} \nonumber \\& \qquad+ \bigg [\left (\delta _{\alpha \delta } + \frac {\boldsymbol{\varLambda }_{\alpha \delta }^2}{a^2\omega ^2R^2}\right )\big (\varPsi _{\delta [\beta ][\gamma ]}[rh_s] + \varPsi _{\delta (\beta )(\gamma )}[h_s/r]\big ) \nonumber \\& \qquad\quad+ \frac {\kappa }{a^2\omega ^2}\big (\varPsi _{|\alpha |[\beta ][\gamma ]}[r] + \varPsi _{|\alpha |(\beta )(\gamma )}[1/r]\big )\bigg ]\sum _b\frac {P_\beta ^{a - b}P_\gamma ^b}{2(a - b)b\omega ^2}\Bigg \}. \end{align}
Combining this with the left-hand side of (2.35a ), the same leftover boundary term vanishes, to give the full projected mass conservation equation
\begin{align} \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s} \right .&\left .{}\!\!- \varPsi _{\{\alpha \}\beta }[r]\right )U_\beta ^a - {\textrm{i}}\kern0.5pt a\omega \left [\left (1 - \frac {\lambda _\alpha ^2}{a^2\omega ^2R^2}\right )\varPsi _{\alpha \beta }[rh_s] - \frac {\kappa }{a^2\omega ^2}\varPsi _{|\alpha |\beta }[r]\right ]P_\beta ^a \nonumber\\=&\, {\textrm{i}}\kern0.5pt a\omega \sum _b\Bigg \{- \beta_0\varPsi _{\alpha \beta \gamma }[rh_s]P_\beta ^{a - b}P_\gamma ^b \nonumber\\&+ \bigg [\left (\delta _{\alpha \delta } + \frac {\boldsymbol{\varLambda }_{\alpha \delta }^2}{a^2\omega ^2R^2}\right )\varPsi _{\delta \beta \gamma }[rh_s] + \frac {\kappa }{a^2\omega ^2}\varPsi _{|\alpha |\beta \gamma }[r]\bigg ]\frac {P_\beta ^{a - b}P_\gamma ^b - U_\beta ^{a - b}U_\gamma ^b}{2} \nonumber\\&+ \bigg [\left (\delta _{\alpha \delta } + \frac {\boldsymbol{\varLambda }_{\alpha \delta }^2}{a^2\omega ^2R^2}\right )\big (\varPsi _{\delta [\beta ][\gamma ]}[rh_s] + \varPsi _{\delta (\beta )(\gamma )}[h_s/r]\big ) \nonumber\\&\quad \quad \quad \quad \quad + \frac {\kappa }{a^2\omega ^2}\big (\varPsi _{|\alpha |[\beta ][\gamma ]}[r] + \varPsi _{|\alpha |(\beta )(\gamma )}[1/r]\big )\bigg ]\dfrac {P_\beta ^{a - b}P_\gamma ^b}{2(a - b)b\omega ^2}\Bigg \}. \end{align}
As before, the left-hand side of the momentum conservation (2.35b ) has no boundary cancellations, so projection is very straightforward
\begin{equation} \iint _S\psi _\alpha \left (\frac {\partial P^a}{\partial s} - {\textrm{i}}\kern0.5pt a\omega h_sU^a\right ){\rm d} S = \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s} + \varPsi _{\alpha \{\beta \}}[r]\right )P_\beta ^a - {\textrm{i}}\kern0.5pt a\omega \varPsi _{\alpha \beta }[rh_s]U_\beta ^a. \end{equation}
Finally, the right-hand side of (2.35b ) requires only that we deal carefully with the tranverse Laplacian, which produces a boundary term
\begin{align} & \iint _S{\textrm{i}}\kern0.5pt\omega h_s\psi _\alpha \sum _bU^{a - b}\left (-b\frac {{\nabla }_{{t}}^2}{b^2\omega ^2}P^b\right ){\rm d} S = {\textrm{i}}\kern0.5pt\omega \int _0^{2 \pi}\left [rh_s\psi _\alpha \sum _b\frac {-b}{b^2\omega ^2}U^{a - b}\frac {\partial P^b}{\partial r}\right ]_0^R{\rm d} \theta \nonumber \\[3pt]& \qquad \qquad- {\textrm{i}}\kern0.5pt\omega \iint _S\sum _b\left (\frac {-b}{b^2\omega ^2}\right )\boldsymbol{\nabla }_{{t}}(h_s\psi _\alpha U^{a - b})\boldsymbol{\cdot }\boldsymbol{\nabla }_{{t}}(\psi _\gamma )P_\gamma ^b{\rm d} S. \end{align}
As took place in two dimensions, the boundary term here becomes a
$U^{a - b}U^b$
, which can be expanded, and the other term has its integration by parts reversed once more, leaving
\begin{align} & \iint _S{\textrm{i}}\kern0.5pt\omega h_s\psi _\alpha \sum _bU^{a - b}\left (-b\frac {{\nabla }_{{t}}^2}{b^2\omega ^2}P^b\right ){\rm d} S = R'\varPsi _{[\alpha \beta \gamma ]}\sum _bU_\beta ^{a - b}U_\gamma ^b \nonumber \\[3pt]& \quad - {\textrm{i}}\kern0.5pt\omega \int _0^{2 \pi}\left [\sum _b\frac {-b}{b^2\omega ^2}rh_s\psi _\alpha U^{a - b}\frac {\partial \psi _\gamma }{\partial r}P_\gamma ^b\right ]_0^R{\rm d} \theta + {\textrm{i}}\kern0.5pt\omega \sum _b\frac {b\boldsymbol{\varLambda }_{\delta \gamma }^2}{b^2\omega ^2R^2}\varPsi _{\alpha \beta \delta }[rh_s] U_\beta ^{a - b}P_\gamma ^b. \end{align}
The leftover boundary term here cancels due to the Neumann condition on the spatial modes, so the full projected version of (2.35b ) is
\begin{align} \left (\delta _{\alpha \beta }\frac {{\rm d} }{{\rm d} s} + \varPsi _{\alpha \{\beta \}}[r]\right )P_\beta ^a - {\textrm{i}}\kern0.5pt a\omega \varPsi _{\alpha \beta }[rh_s]U_\beta ^a = R'\varPsi _{[\alpha \beta \gamma ]}\sum _bU_\beta ^{a - b}U_\gamma ^b \nonumber \\ + {\textrm{i}}\kern0.5pt\omega \sum _b\Bigg \{\varPsi _{\alpha \beta \delta }[rh_s]\left [(a - b)\delta _{\delta \gamma } - b\left (\delta _{\delta \gamma } - \frac {\boldsymbol{\varLambda }_{\delta \gamma }^2}{b^2\omega ^2R^2}\right )\right ] \nonumber \\ + \frac {\varPsi _{\alpha [\beta ][\gamma ]}[rh_s] + \varPsi _{\alpha (\beta )(\gamma )}[h_s/r]}{b\omega ^2}\Bigg \}U_\beta ^{a - b}P_\gamma ^b. \end{align}
We now have the same task of rewriting these equations in a more numerically efficient way. The added complication in three dimensions is that these are double integrals, which means that they must be split into separate radial and azimuthal parts; nonetheless, everything may still be made independent of
$s$
just as before. The notation for the split integrals is very similar to that of the
$\varPsi$
matrices, e.g.
\begin{equation} \varPi _{\alpha [\beta ]}[f(x)] = C_\alpha C_\beta \int _{x = 0}^1{{J}}_{p_\alpha }(\lambda _\alpha x)\frac {\partial }{\partial x}\left [{{J}}_{p_\beta }(\lambda _\beta x)\right ] f(x)\,{\rm d} x, \end{equation}
and
\begin{equation} \varPhi _{\alpha (\beta )}[g(\phi )] = \frac {1}{ \pi}\int _{\theta = 0}^{2 \pi}\cos \left(p_\alpha \phi - \xi _\alpha \frac { \pi}{2} \right)\frac {\partial }{\partial \theta }\left [\cos \left(p_\beta \phi - \xi _\beta \frac { \pi}{2} \right)\right ]g(\phi )\,{\rm d} \theta . \end{equation}
Once again, we need to know how
$s$
-derivatives of modes work: in three dimensions, the relevant equation is
\begin{equation} \frac {\partial \psi _\alpha }{\partial s} = -\frac {R'}{R}\left (\psi _\alpha + r\frac {\partial \psi _\alpha }{\partial r}\right ) - \tau \frac {\partial \psi _\alpha }{\partial \theta }, \end{equation}
while the analogous improved forms of matrices with derivatives on multiple subscripts are
\begin{align} & \varPsi _{\alpha [\beta ][\gamma ]}[rh_s] + \varPsi _{\alpha (\beta )(\gamma )}[h_s/r] \nonumber \\[3pt]& \qquad = \frac {\lambda _\beta ^2 + \lambda _\gamma ^2 - \lambda _\alpha ^2}{2R^2}\varPsi _{\alpha \beta \gamma }[rh_s] - \kappa \varPsi _{|\alpha |\beta \gamma }[r] + \frac {\kappa }{2}\varPsi _{[\alpha \beta \gamma ]}[r\cos \phi ], \\[-12pt] \nonumber \end{align}
\begin{align} & \varPsi _{|\alpha |[\beta ][\gamma ]}[r] + \varPsi _{|\alpha |(\beta )(\gamma )}[1/r] \nonumber \\[3pt]& \qquad = \frac {\lambda _\beta ^2 + \lambda _\gamma ^2 - \lambda _\alpha ^2}{2R^2}\varPsi _{|\alpha |\beta \gamma }[r] + \frac {\lambda _\alpha ^2 - p_\alpha ^2}{2R^2}\varPsi _{[\alpha \beta \gamma ]}[r\cos \phi ] - \frac {1}{2R}\varPsi _{[(\alpha )\beta \gamma ]}[\sin \phi ]. \\[9pt] \nonumber \end{align}
Equipped with the above, (A25) and (A29) become
\begin{align}& \left (\dfrac {{\rm d} }{{\rm d} s} + \dfrac {R'}{R} \unicode{x1D652} + \tau \unicode{x1D643}\right )\boldsymbol{u}^a - {\textrm{i}}\kern0.5pt a\omega \left [\left ( \unicode{x1D644} - \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R^2}\right )\left ( \unicode{x1D644} - \kappa R \unicode{x1D63C}\right ) - \dfrac {\kappa \kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}}{a^2\omega ^2R}\right ]\boldsymbol{p}^a \nonumber \\& = \dfrac {{\textrm{i}}\kern0.5pt a\omega }{\sqrt { \pi}R}\sum _b\left \{\vphantom{\left [\dfrac {\left ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R^2}\right )\left (\mathcal{I}^\lambda - \kappa R\mathcal{A}^\lambda \right ) + \dfrac {\kappa \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda }{a^2\omega ^2R}}{2(a - b)b\omega ^2R^2} - \beta_0\left (\mathcal{I} - \kappa R\mathcal{A}\right )\right ]}\left [\left ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R^2}\right )\left (\mathcal{I} - \kappa R\mathcal{A}\right ) + \dfrac {\kappa \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}}{a^2\omega ^2 R}\right ]\dfrac {\langle \boldsymbol{p}^{a - b},\boldsymbol{p}^b\rangle - \langle \boldsymbol{u}^{a - b},\boldsymbol{u}^b\rangle }{2} \right.\nonumber \\& \left.\qquad + \left [\dfrac {\left ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R^2}\right )\left (\mathcal{I}^\lambda - \kappa R\mathcal{A}^\lambda \right ) + \dfrac {\kappa \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda }{a^2\omega ^2R}}{2(a - b)b\omega ^2R^2} - \beta_0\left (\mathcal{I} - \kappa R\mathcal{A}\right )\right ]\langle \boldsymbol{p}^{a - b},\boldsymbol{p}^b\rangle \right \}, \\[-12pt] \nonumber\end{align}
\begin{align} & \left (\frac {{\rm d} }{{\rm d} s} - \frac {R'}{R} \unicode{x1D652}^{{\,T}} - \tau \unicode{x1D643}^{{T}}\right )\boldsymbol{p}^a - {\textrm{i}}\kern0.5pt a\omega \big ( \unicode{x1D644} - \kappa R \unicode{x1D63C}\big )\boldsymbol{u}^a = \frac {1}{\sqrt { \pi}R}\sum _b\left \{\vphantom{\left [\left (\mathcal{I} - \kappa R\mathcal{A}\right )\bigg \langle \unicode{x1D644},(a - b) \unicode{x1D644} - b\left ( \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{b^2\omega ^2R^2}\right )\bigg \rangle + \frac {\mathcal{I}^\lambda - \kappa R\mathcal{A}^\lambda }{b\omega ^2R^2}\right ]}\frac {R'}{R}\overline {\mathcal{W}}\langle \boldsymbol{u}^{a - b},\boldsymbol{u}^b\rangle \nonumber\right. \\& \left.\quad + {\textrm{i}}\kern0.5pt\omega \left [\left (\mathcal{I} - \kappa R\mathcal{A}\right )\bigg \langle \unicode{x1D644},(a - b) \unicode{x1D644} - b\left ( \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{b^2\omega ^2R^2}\right )\bigg \rangle + \frac {\mathcal{I}^\lambda - \kappa R\mathcal{A}^\lambda }{b\omega ^2R^2}\right ]\langle \boldsymbol{u}^{a - b},\boldsymbol{p}^b\rangle \right \}, \\[9pt] \nonumber \end{align}
where
\begin{align} \boldsymbol{\varLambda }_{\alpha \beta } &= \lambda _\alpha \delta _{\alpha \beta }, & & \unicode{x1D64B}_{\alpha \beta } = p_\alpha \delta _{\alpha \beta }, \\[-12pt] \nonumber \end{align}
\begin{align} \unicode{x1D652}_{\alpha \beta } &= \delta _{\alpha \beta } + \varPi _{[\alpha ]\beta }[x^2]\varPhi _{\alpha \beta }, & & \unicode{x1D643}_{\alpha \beta } = \varPi _{\alpha \beta }[x]\varPhi _{(\alpha )\beta }, \\[-12pt] \nonumber \end{align}
\begin{align} \unicode{x1D63C}_{\alpha \beta } &= \varPi _{\alpha \beta }[x^2]\varPhi _{\alpha \beta }[\cos \phi ], & &\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}_{\alpha \beta } = \varPi _{[\alpha ]\beta }[x]\varPhi _{\alpha \beta }[\cos \phi ] - \varPi _{\alpha \beta }\varPhi _{(\alpha )\beta }[\sin \phi ], \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{I}_{\alpha \beta \gamma } &= \varPi _{\alpha \beta \gamma }[x]\varPhi _{\alpha \beta \gamma }, & &\overline {\mathcal{W}}_{\alpha \beta \gamma } = \varPi _{[\alpha \beta \gamma ]}\varPhi _{\alpha \beta \gamma }, \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{A}_{\alpha \beta \gamma } &= \varPi _{\alpha \beta \gamma }[x^2]\varPhi _{\alpha \beta \gamma }[\cos \phi ], & &\kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}_{\alpha \beta \gamma } = \varPi _{[\alpha ]\beta \gamma }[x]\varPhi _{\alpha \beta \gamma }[\cos \phi ] - \varPi _{\alpha \beta \gamma }\varPhi _{(\alpha )\beta \gamma }[\sin \phi ], \\[-12pt] \nonumber \end{align}
\begin{align} \kern2.5pt \overline {\kern-2.5pt \mathcal{A}}_{\alpha \beta \gamma } &= \varPi _{[\alpha \beta \gamma ]}\varPhi _{\alpha \beta \gamma }[\cos \phi ], & &\kern2.5pt \overline {\kern-2.5pt \mathcal{A}}_{\alpha \beta \gamma }^* = \varPi _{[\alpha \beta \gamma ]}\varPhi _{(\alpha )\beta \gamma }[\sin \phi ], \\[9pt] \nonumber \end{align}
and the eigenvalue-modified matrices
$\mathcal{I}^\lambda$
,
$\mathcal{A}^\lambda$
and
$\kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda$
are defined by
\begin{align} \mathcal{I}^\lambda &= \mathcal{I}\frac {\langle \boldsymbol{\varLambda }^2, \unicode{x1D644}\,\rangle + \langle \unicode{x1D644},\boldsymbol{\varLambda }^2\rangle - \boldsymbol{\varLambda }^2}{2}, \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{A}^\lambda &= \mathcal{A}\frac {\langle \boldsymbol{\varLambda }^2, \unicode{x1D644}\,\rangle + \langle \unicode{x1D644},\boldsymbol{\varLambda }^2\rangle - \boldsymbol{\varLambda }^2}{2} + \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}} - \frac {\kern2.5pt \overline {\kern-2.5pt \mathcal{A}}}{2}, \\[-12pt] \nonumber \end{align}
\begin{align} \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda &= \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}\frac {\langle \boldsymbol{\varLambda }^2, \unicode{x1D644}\,\rangle + \langle \unicode{x1D644},\boldsymbol{\varLambda }^2\rangle - \boldsymbol{\varLambda }^2}{2} + \frac {(\boldsymbol{\varLambda }^2 - \unicode{x1D64B}^{\,2})\kern2.5pt \overline {\kern-2.5pt \mathcal{A}} - \kern2.5pt \overline {\kern-2.5pt \mathcal{A}}^*}{2}. \\[9pt] \nonumber \end{align}
As before, notation has been chosen such that curvature acts through annularity matrices (
$ \unicode{x1D63C}$
etc.), width variation through
$ \unicode{x1D652}$
, and with the new quantity, torsion, through a helicity matrix
$ \unicode{x1D643}$
. (A34) are (2.49) in the main text.
A.3. Evaluation of integrals
Here, the mode integrals employed in the projection of the governing equations above are calculated explicitly, where possible. These are the expressions used when performing computations, in order to save time.
A.3.1. Integrals from the two-dimensional formulation
We have
\begin{align} \varXi _{\alpha \beta } &= \delta _{\alpha \beta }, \\[-12pt] \nonumber \end{align}
\begin{align} \varXi _{\alpha \beta }[\xi ] &= \frac {\delta _{\alpha \beta }}{2} + C_\alpha C_\beta \frac {\left ((-1)^{\alpha + \beta } - 1\right )(\alpha ^2 + \beta ^2)}{(\alpha ^2 - \beta ^2)^2 \pi^2}, \\[-12pt] \nonumber \end{align}
\begin{align} \varXi _{[\alpha ]\beta } &= \sqrt {2}\alpha ^2 C_\beta \frac {(-1)^{\alpha + \beta } - 1}{\alpha ^2 - \beta ^2}, \\[-12pt] \nonumber \end{align}
\begin{align} \varXi _{[\alpha ]\beta }[\xi ] &= \delta _{\alpha \beta }\frac {(1 - \delta _{\alpha 0})}{2} + (1 - \delta _{\alpha \beta })\frac {\sqrt {2}\alpha ^2C_\beta (-1)^{\alpha + \beta }}{\alpha ^2 - \beta ^2}, \\[-12pt] \nonumber \end{align}
\begin{align} \varXi _{\alpha \beta \gamma } &= \frac {C_\alpha C_\beta }{2 C_\gamma }(\delta _{\alpha +\beta ,\gamma } + \delta _{|\alpha - \beta |,\gamma }), \\[-12pt] \nonumber \end{align}
\begin{align} \varXi _{\alpha \beta \gamma }[\xi ] &= \frac {\varXi _{\alpha \beta \gamma }}{2} + \frac {C_\alpha C_\beta C_\gamma \left ((-1)^{\alpha + \beta + \gamma } - 1\right )}{4 \pi^2}\Bigg (\frac {1}{(\alpha + \beta + \gamma )^2}\notag \\ &\quad \quad \quad + \frac {1}{(\alpha + \beta - \gamma )^2} + \frac {1}{(\alpha - \beta + \gamma )^2} + \frac {1}{(\alpha - \beta + \gamma )^2}\Bigg ), \\[-12pt] \nonumber \end{align}
\begin{align} \varXi _{[\alpha ]\beta \gamma } &= \frac {\alpha C_\beta C_\gamma \left ((-1)^{\alpha + \beta + \gamma } - 1\right )}{2\sqrt {2}}\Bigg (\frac {1}{\alpha + \beta - \gamma } + \frac {1}{\alpha - \beta + \gamma } + \frac {1}{\alpha - \beta - \gamma }\Bigg ), \\[-12pt] \nonumber \end{align}
\begin{align} \varXi _{[\alpha \beta \gamma ]} &= C_\alpha C_\beta C_\gamma \left ((-1)^{\alpha + \beta + \gamma } - 1\right ), \\[-12pt] \nonumber \end{align}
\begin{align} \varXi _{[\alpha \beta \gamma ]}^+ &= C_\alpha C_\beta C_\gamma (-1)^{\alpha + \beta + \gamma }, \\[-12pt] \nonumber \end{align}
\begin{align} \varXi _{[\alpha \beta \gamma ]}^- &= C_\alpha C_\beta C_\gamma . \\[9pt] \nonumber \end{align}
A.3.2. Integrals from the three-dimensional formulation
We have
\begin{align} \varPi _{\alpha \beta } &= C_\alpha C_\beta \int _0^1 J_{p_\alpha }(\lambda _\alpha x)J_{p_\beta }(\lambda _\beta x){\rm d}x, \\[-12pt] \nonumber \end{align}
\begin{align} \varPi _{\alpha \beta }[x] &= C_\alpha C_\beta \int _0^1xJ_{p_\alpha }(\lambda _\alpha x)J_{p_\beta }(\lambda _\beta x){\rm d}x, \\[-12pt] \nonumber \end{align}
\begin{align} \varPi _{\alpha \beta }[x^2] &= C_\alpha C_\beta \int _0^1 x^2J_{p_\alpha }(\lambda _\alpha x)J_{p_\beta }(\lambda _\beta x){\rm d}x, \\[-12pt] \nonumber \end{align}
\begin{align} \varPi _{[\alpha ]\beta }[x] &= C_\alpha C_\beta \int _0^1 \left [p_\alpha J_{p_\alpha }(\lambda _\alpha x) - \lambda _\alpha xJ_{p_\alpha + 1}(\lambda _\alpha x)\right ]J_{p_\beta }(\lambda _\beta x){\rm d}x, \\[-12pt] \nonumber \end{align}
\begin{align} \varPi _{[\alpha ]\beta }[x^2] &= C_\alpha C_\beta \int _0^1 x\left [p_\alpha J_{p_\alpha }(\lambda _\alpha x) - \lambda _\alpha xJ_{p_\alpha + 1}(\lambda _\alpha x)\right ]J_{p_\beta }(\lambda _\beta x){\rm d}x, \\[-12pt] \nonumber \end{align}
\begin{align} \varPhi _{\alpha \beta } &= (\delta _{p_\alpha + p_\beta ,0} + \delta _{p_\alpha p_\beta })\delta _{\xi _\alpha \xi _\beta }, \\[-12pt] \nonumber \end{align}
\begin{align} \varPhi _{(\alpha )\beta } &= -p_\alpha (-1)^{\xi _\alpha }\delta _{p_\alpha p_\beta }\delta _{\xi _\alpha + \xi _\beta ,1}, \\[-12pt] \nonumber \end{align}
\begin{align} \varPhi _{\alpha \beta }[\cos \phi ] &= \frac {1}{2}\left [(-1)^{\xi _\alpha }\delta _{p_\alpha + p_\beta ,1} + \delta _{|p_\alpha - p_\beta |,1}\right ]\delta _{\xi _\alpha \xi _\beta }, \\[-12pt] \nonumber \end{align}
\begin{align} \varPhi _{\alpha \beta }[\sin \phi ] &= \frac {1}{2}\left [\delta _{p_\alpha + p_\beta ,1} + (-1)^{\xi _\beta }(p_\alpha - p_\beta )\delta _{|p_\alpha - p_\beta |,1}\right ]\delta _{\xi _\alpha + \xi _\beta ,1}, \\[-12pt] \nonumber \end{align}
\begin{align} \varPhi _{(\alpha )\beta }[\sin \phi ] &= -\frac {p_\alpha }{2}\left [\delta _{p_\alpha 1}\delta _{p_\beta 0}\delta _{\xi _\alpha 0} + (p_\alpha - p_\beta )\delta _{|p_\alpha - p_\beta |,1}\right ]\delta _{\xi _\alpha \xi _\beta }. \\[9pt] \nonumber \end{align}
The linear coefficient matrices in three dimensions always consist of a combination of a
$\varPi$
with a
$\varPhi$
, which rarely simplifies further, although the notable cases are the identity
\begin{equation} \mathsf{I}_{\alpha \beta } = \varPi _{\alpha \beta }[x]\varPhi _{\alpha \beta } = \delta _{\alpha \beta }, \end{equation}
and the helicity matrix
\begin{equation} \mathsf{H}_{\alpha \beta } = \varPi _{\alpha \beta }[x]\varPhi _{(\alpha )\beta } = -p_\alpha (-1)^{\xi _\alpha }\delta _{\xi _\alpha +\xi _\beta ,1}\delta _{p_\alpha p_\beta }\delta _{q_\alpha q_\beta }. \end{equation}
Likewise,
\begin{align} \varPi _{\alpha \beta \gamma } &= C_\alpha C_\beta C_\gamma \!\!\int _0^1 J_{p_\alpha }(\lambda _\alpha x)J_{p_\beta }(\lambda _\beta x)J_{p_\gamma }(\lambda _\gamma x){\rm d}x, \\[-12pt] \nonumber \end{align}
\begin{align} \varPi _{\alpha \beta \gamma }[x] &= C_\alpha C_\beta C_\gamma \!\!\int _0^1 xJ_{p_\alpha }(\lambda _\alpha x)J_{p_\beta }(\lambda _\beta x)J_{p_\gamma }(\lambda _\gamma x){\rm d}x, \\[-12pt] \nonumber \end{align}
\begin{align} \varPi _{\alpha \beta \gamma }[x^2] &= C_\alpha C_\beta C_\gamma \!\!\int _0^1 x^2J_{p_\alpha }(\lambda _\alpha x)J_{p_\beta }(\lambda _\beta x)J_{p_\gamma }(\lambda _\gamma x){\rm d}x, \\[-12pt] \nonumber \end{align}
\begin{align} \varPi _{[\alpha ]\beta \gamma }[x] &= C_\alpha C_\beta C_\gamma \!\!\int _0^1 \big (p_\alpha J_{p_\alpha }(\lambda _\alpha x) - \lambda _\alpha xJ_{p_\alpha + 1}(\lambda _\alpha x)\big )J_{p_\beta }(\lambda _\beta x)J_{p_\gamma }(\lambda _\gamma x){\rm d}x, \\[-12pt] \nonumber\end{align}
\begin{align} \varPi _{[\alpha \beta \gamma ]}[x] &= C_\alpha C_\beta C_\gamma J_{p_\alpha }(\lambda _\alpha )J_{p_\beta }(\lambda _\beta )J_{p_\gamma }(\lambda _\gamma ), \\[-12pt] \nonumber \end{align}
\begin{align} \varPhi _{\alpha \beta \gamma } &= \frac {1}{2}\big (\delta _{p_\alpha + p_\beta + p_\gamma ,0} + (-1)^{\xi _\alpha }\delta _{p_\beta + p_\gamma ,p_\alpha } + (-1)^{\xi _\beta }\delta _{p_\alpha + p_\gamma ,p_\beta }\notag \\ & \quad + (-1)^{\xi _\gamma }\delta _{p_\alpha + p_\beta ,p_\gamma }\big )\big (\delta _{\xi _\alpha + \xi _\beta + \xi _\gamma ,0} - \delta _{\xi _\alpha + \xi _\beta + \xi _\gamma ,2}\big ), \\[-12pt] \nonumber \end{align}
\begin{align} \varPhi _{\alpha \beta \gamma }[\cos \phi ] &= \frac {1}{4}\big (\delta _{p_\alpha + p_\beta + p_\gamma ,1} + (-1)^{\xi _\alpha }\delta _{|p_\beta + p_\gamma - p_\alpha |,1} + (-1)^{\xi _\beta }\delta _{|p_\alpha + p_\gamma - p_\beta |,1}\notag \\ &\quad + (-1)^{\xi _\gamma }\delta _{|p_\alpha + p_\beta - p_\gamma |,1}\big )\big (\delta _{\xi _\alpha + \xi _\beta + \xi _\gamma ,0} - \delta _{\xi _\alpha + \xi _\beta + \xi _\gamma ,2}\big ), \\[-12pt] \nonumber \end{align}
\begin{align} \varPhi _{(\alpha )\beta \gamma }[\sin \phi ] &= -\frac {p_\alpha }{4}\big (\delta _{p_\alpha + p_\beta + p_\gamma ,1} + \delta _{|p_\beta + p_\gamma - p_\alpha |,1}(-1)^{p_\beta + p_\gamma - p_\alpha + \xi _\alpha }\notag \\ &\quad + \delta _{|p_\alpha + p_\gamma - p_\beta |,1}(-1)^{p_\alpha + p_\gamma - p_\beta + \xi _\beta } + \delta _{|p_\alpha + p_\beta - p_\gamma |,1}(-1)^{p_\alpha + p_\beta - p_\gamma + \xi _\gamma }\big ]\big )\notag \\ &\quad\times \big (\delta _{\xi _\alpha + \xi _\beta + \xi _\gamma ,0} - \delta _{\xi _\alpha + \xi _\beta + \xi _\gamma ,2}\big ). \\[9pt] \nonumber \end{align}
Appendix B. Definitions of block matrices
In this appendix we detail the entries of the block matrix
$ \unicode{x1D647}^a$
and block tensor
$\mathcal{N}^{\,ab}$
from (2.50). In terms of square blocks
$ \unicode{x1D647}^a =: \begin{pmatrix} \unicode{x1D647}_1^a & \unicode{x1D647}_2^a \\[5pt] \unicode{x1D647}_3^a & \unicode{x1D647}_4^a \end{pmatrix}$
, we have a two-dimensional
$ \unicode{x1D647}^a$
\begin{equation} \unicode{x1D647}^a = \begin{pmatrix} -\dfrac {X'}{2X} \unicode{x1D652} - \dfrac {X_-'}{X}\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}} &{\textrm{i}}\kern0.5pt a\omega \left [\left ( \unicode{x1D644} - \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X^2}\right )\left ( \unicode{x1D644} - \kappa X \unicode{x1D63C}\right ) - \dfrac {\kappa \kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}}{a^2\omega ^2X}\right ] \\[9pt] {\textrm{i}}\kern0.5pt a\omega ( \unicode{x1D644} - \kappa X \unicode{x1D63C}) &\dfrac {X'}{2X} \unicode{x1D652}^T + \dfrac {X_-'}{X}\kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}^T \end{pmatrix}\!, \end{equation}
and a three-dimensional
$ \unicode{x1D647}^a$
\begin{equation} \unicode{x1D647}^a = \begin{pmatrix} -\dfrac {R'}{R} \unicode{x1D652} - \tau \unicode{x1D643} &{\textrm{i}}\kern0.5pt a\omega \left [\left ( \unicode{x1D644} - \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R^2}\right )\left ( \unicode{x1D644} - \kappa R \unicode{x1D63C}\right ) - \dfrac {\kappa \kern1.5pt\widetilde {\kern-1.5pt \unicode{x1D63C}}}{a^2\omega ^2R}\right ] \\[5pt] {\textrm{i}}\kern0.5pt a\omega ( \unicode{x1D644} - \kappa R \unicode{x1D63C}) &\dfrac {R'}{R} \unicode{x1D652}^T + \tau \unicode{x1D643}^T \end{pmatrix}\!, \end{equation}
and in terms of ‘cubic blocks’
\begin{align} &\mathcal{N}_{\,\alpha ,\beta ,0:\alpha _{\textrm{max}}}^{\,ab} = \begin{pmatrix} \mathcal{N}_1^{ab} &\mathcal{N}_2^{ab} \\[5pt] \mathcal{N}_3^{ab} &\mathcal{N}_4^{ab} \end{pmatrix}, && \mathcal{N}_{\,\alpha ,\beta ,\alpha _{\textrm{max}}+1:2\alpha _{\textrm{max}}+1}^{\,ab} = \begin{pmatrix} \mathcal{N}_5^{ab} &\mathcal{N}_6^{ab} \\[5pt] \mathcal{N}_7^{ab} &\mathcal{N}_8^{ab} \end{pmatrix}, \end{align}
we have, in two dimensions
\begin{align} \mathcal{N}_1^{ab} &= -\frac {{\textrm{i}}\kern0.5pt a\omega }{2\sqrt {X}}\left [\left ( \unicode{x1D644} + \frac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X^2}\right )\left (\mathcal{I} - \kappa X\mathcal{A}\right ) + \frac {\kappa \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}}{a^2\omega ^2 X}\right ], \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{N}_3^{ab} &= \frac {1}{\sqrt {X}}\bigg (\frac {X_+'}{X}\overline {\mathcal{W}}^+ - \frac {X_-'}{X}\overline {\mathcal{W}}^-\bigg ), \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{N}_6^{ab} &= -\mathcal{N}_1^{ab} + \dfrac {{\textrm{i}}\kern0.5pt a\omega }{\sqrt {X}}\left [\dfrac {\left ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2X^2}\right )\left (\mathcal{I}^\lambda - \kappa X\mathcal{A}^\lambda \right ) + \dfrac {\kappa \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda }{a^2\omega ^2X}}{2(a - b)b\omega ^2X^2} - \beta_0\left (\mathcal{I} - \kappa R\mathcal{A}\right )\right ]\!, \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{N}_7^{ab} &= \frac {{{i}} \omega }{\sqrt {X}}\left \{\left (\mathcal{I} - \kappa X\mathcal{A}\right )\bigg \langle \unicode{x1D644},(a - b) \unicode{x1D644} - b\left ( \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{b^2\omega ^2X^2}\right )\bigg \rangle + \frac {\mathcal{I}^\lambda - \kappa X\mathcal{A}^\lambda }{b\omega ^2X^2}\right \}\!, \\[9pt] \nonumber \end{align}
and in three dimensions
\begin{align} \mathcal{N}_1^{ab} &= -\frac {{\textrm{i}}\kern0.5pt a\omega }{2\sqrt { \pi}R}\left [\left ( \unicode{x1D644} + \frac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R^2}\right )\left (\mathcal{I} - \kappa R\mathcal{A}\right ) + \frac {\kappa \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}}{a^2\omega ^2 R}\right ]\!, \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{N}_3^{ab} &= \frac {R'}{\sqrt { \pi}R^2}\overline {\mathcal{W}}, \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{N}_6^{ab} &= -\mathcal{N}_1^{ab} + \dfrac {{\textrm{i}}\kern0.5pt a\omega }{\sqrt { \pi}R}\left [\dfrac {\left ( \unicode{x1D644} + \dfrac {\boldsymbol{\varLambda }^2}{a^2\omega ^2R^2}\right )\left (\mathcal{I}^\lambda - \kappa R\mathcal{A}^\lambda \right ) + \dfrac {\kappa \kern2.5pt \widetilde {\kern-2.5pt \mathcal{A}}^\lambda }{a^2\omega ^2R}}{2(a - b)b\omega ^2R^2} - \beta_0\left (\mathcal{I} - \kappa R\mathcal{A}\right )\right ]\!, \\[-12pt] \nonumber \end{align}
\begin{align} \mathcal{N}_7^{ab} &= \frac {{{i}} \omega }{\sqrt { \pi}R}\left \{\left (\mathcal{I} - \kappa R\mathcal{A}\right )\bigg \langle \unicode{x1D644},(a - b) \unicode{x1D644} - b\left ( \unicode{x1D644} - \frac {\boldsymbol{\varLambda }^2}{b^2\omega ^2R^2}\right )\bigg \rangle + \frac {\mathcal{I}^\lambda - \kappa R\mathcal{A}^\lambda }{b\omega ^2R^2}\right \}\!, \\[9pt] \nonumber \end{align}
while in both two and three dimensions
${\mathcal{N}}_2^{ab} = {\mathcal{N}}_4^{ab} = {\mathcal{N}}_5^{ab} = {\mathcal{N}}_8^{ab} = 0$
.
Appendix C. Characteristic admittances of non-straight ducts
This appendix deals with the non-straight-duct versions of § 3.2.1.
C.1. Curved-duct characteristic admittances
When we introduce curvature, the straight-duct operator from (3.5) retains its block structure
\begin{equation} \breve { \unicode{x1D647}}^a = \unicode{x1D647}^a\Bigg |_{\textrm{no}\,\tau ,\,X'\,{or}\,R'} = \begin{pmatrix} \boldsymbol{\mathit{0}} &\breve { \unicode{x1D647}}_2^a \\[4pt] \breve { \unicode{x1D647}}_{\,3}^{\,a} &\boldsymbol{\mathit{0}} \end{pmatrix}\!, \end{equation}
with the difference that
$\breve { \unicode{x1D647}}_2^a$
and
$\breve { \unicode{x1D647}}_{\,3}^{\,a}$
are no longer diagonal. They do remain symmetric, however (by inspection in
$\breve { \unicode{x1D647}}_{\,3}^{\,a}$
’s case, and by a complicated integration-by-parts in
$\breve { \unicode{x1D647}}_2^a$
’s). The eigenvalues
$\breve {\gamma }^a$
now satisfy (from (3.7))
\begin{equation} 0 = \det \left [\breve { \unicode{x1D647}}^a - \breve {\gamma }^a\begin{pmatrix} \unicode{x1D644} &\boldsymbol{\mathit{0}} \\ \boldsymbol{\mathit{0}} & \unicode{x1D644} \end{pmatrix}\right ] = \det \left ((\breve {\gamma }^{a})^2 \unicode{x1D644} - \breve { \unicode{x1D647}}_{\,3}^{\,a}\breve { \unicode{x1D647}}_2^a\right )\!, \end{equation}
so they will exhibit the same properties as before (purely real/imaginary and being mirrored either side of the origin). We once again have
$\breve { \unicode{x1D647}}_{\,3}^{\,a}\breve { \unicode{x1D647}}_2^a\breve {\boldsymbol{c}}_{\alpha ,2}^{a\pm } = (\breve {\gamma }_{\,\alpha}^{\,a})^2\breve {\boldsymbol{c}}_{\alpha ,2}^{a\pm }$
, so a single lower-eigenvector matrix
$\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^a$
still exists, but this time it is non-diagonal. Substituting the admittance into the lower block this time, we get
\begin{align} &\breve { \unicode{x1D647}}_{\,3}^{\,a}\kern2.5pt \breve { \kern-2.5pt \unicode{x1D654}}^{\,a\pm }\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^a = \pm \kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^a\breve {\boldsymbol{\varGamma }}^a,&& \textrm{so}&& \kern2.5pt \breve { \kern-2.5pt \unicode{x1D654}}^{\,a\pm } = \pm \left(\breve { \unicode{x1D647}}_{\,3}^{\,a} \right)^{-1}\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^a\breve {\boldsymbol{\varGamma }}^a \big(\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^a \big)^{-1} =:\, \pm \breve { \kern-3pt \unicode{x1D654}}^{\,a}. \end{align}
The nonlinear characteristic admittance takes a bit more work this time: the term multiplying
$\breve {\mathcal{Y}}^{\,ab\pm }$
in the nonlinear admittance equation is no longer diagonal. We have already diagonalised
$\breve { \unicode{x1D647}}_{\,3}^{\,a}\kern2.5pt \breve { \kern-2.5pt \unicode{x1D654}}^{\,a\pm }$
, but we also need to diagonalise
$\kern2.5pt \breve { \kern-2.5pt \unicode{x1D654}}^{\,a\pm }\breve { \unicode{x1D647}}_{\,3}^{\,a}$
. To do this we note that (3.2a) is satisfied by
$( \unicode{x1D654}^{\,a})^{{T}}$
as well as
$ \unicode{x1D654}^{\,a}$
, so with symmetric
$\breve { \unicode{x1D647}}_{\,3}^{\,a}$
, we have
\begin{equation} \kern2.5pt \breve { \kern-2.5pt \unicode{x1D654}}^{\,a\pm }\breve { \unicode{x1D647}}_{\,3}^{\,a} = (\breve { \unicode{x1D647}}_{\,3}^{\,a}\kern2.5pt \breve { \kern-2.5pt \unicode{x1D654}}^{\,a\pm })^{\textrm{T}} = \pm \big (\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^a\big )^{-{\textrm{T}}}\breve {\boldsymbol{\varGamma }}^a\big (\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^a\big )^{\textrm{T}}. \end{equation}
With both diagonalisations, we now have
\begin{align} \breve {\mathcal{Y}}_{\,\alpha \beta \gamma }^{\,ab\pm } & = \pm \sum _{\zeta ,\delta ,\epsilon }\big (\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^a\big )_{\zeta \alpha }^{-1}\big (\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^{a - b}\big )_{\delta \beta }^{-1}\big (\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^b\big )_{\epsilon \gamma }^{-1} \nonumber \\ &\quad\times \frac {\left [\big (\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^a\big )^{{{T}}}\left [\breve {\mathcal{N}}_1^{\,\,ab}\big\langle \kern2.5pt \breve { \kern-2.5pt \unicode{x1D654}}^{\,\,a - b},\kern2.5pt \breve { \kern-2.5pt \unicode{x1D654}}^b \big\rangle + \breve {\mathcal{N}}_6^{\,\,ab} - \breve { \kern-3pt \unicode{x1D654}}^{\,a}\breve {\mathcal{N}}_{\,\,7}^{\,\,ab} \big\langle \kern2.5pt \breve { \kern-2.5pt \unicode{x1D654}}^{\,\,a - b}, \unicode{x1D644}\kern2pt \big\rangle \right ]\big \langle \kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^{a - b},\kern2.5pt \breve {\kern-2.5pt \unicode{x1D63E}}^b\big \rangle \right ]_{\zeta \delta \epsilon }}{\breve {\gamma }_\zeta ^a + \breve {\gamma }_\delta ^{a - b} + \breve {\gamma }_\epsilon ^b}. \end{align}
C.2. Torsional-duct characteristic admittances
The final parameter that we can switch on is torsion, and this is only possible in the three-dimensional case. When we do so, the matrix
$ \unicode{x1D647}^a$
becomes
\begin{equation} \widetilde { \unicode{x1D647}}^a = \unicode{x1D647}^a\bigg |_{\textrm{no}\,R'} = \begin{pmatrix} \widetilde { \unicode{x1D647}}_1^a &\widetilde { \unicode{x1D647}}_2^a \\[6pt] \widetilde { \unicode{x1D647}}_{\,3}^{\,a} &\widetilde { \unicode{x1D647}}_1^a \end{pmatrix}, \end{equation}
where
$\widetilde { \unicode{x1D647}}_2^a$
and
$\widetilde { \unicode{x1D647}}_{\,3}^{\,a}$
remain real and symmetric, and
$\widetilde { \unicode{x1D647}}_1^a$
is a real antisymmetric matrix (apparent by direct calculation of the matrix
$ \unicode{x1D643}$
). Since a linear combination of the identity and an antisymmetric matrix is non-singular, we know that the characteristic equation reads
\begin{equation} 0 = \det \left (\big [\widetilde { \unicode{x1D647}}_1^a - \widetilde {\gamma }^{a} \unicode{x1D644}\big ]^2 - \big [\widetilde { \unicode{x1D647}}_1^a - \widetilde {\gamma }^{a} \unicode{x1D644}\big ]\widetilde { \unicode{x1D647}}_{\,3}^{\,a}\big [\widetilde { \unicode{x1D647}}_1^a - \widetilde {\gamma }^{a} \unicode{x1D644}\big ]^{-1}\widetilde { \unicode{x1D647}}_2^a\right ), \end{equation}
but also, from the alternative expression in (3.7), that we have
\begin{equation} 0 = \det \left (\big [\widetilde { \unicode{x1D647}}_1^a - \widetilde {\gamma }^{a} \unicode{x1D644}\big ]^2 - \widetilde { \unicode{x1D647}}_2^a\big [\widetilde { \unicode{x1D647}}_1^a - \widetilde {\gamma }^{a} \unicode{x1D644}\big ]^{-1}\widetilde { \unicode{x1D647}}_{\,3}^{\,a}\big [\widetilde { \unicode{x1D647}}_1^a - \widetilde {\gamma }^{a} \unicode{x1D644}\big ]\right ). \end{equation}
The expression on the right will have the same determinant if we take its transpose. Doing so reverses the signs on every term except those with a
$\widetilde {\gamma }^a \unicode{x1D644}$
, demonstrating that the eigenvalue mirroring property remains, even if the eigenvalues are no longer necessarily on the real/imaginary axes. We may still partition them into two matched groups, so we define
\begin{equation} \{\widetilde {\gamma }_{\,\alpha}^{\,a}\}_{\alpha = 0}^\infty = \{\widetilde {\gamma }^a : \textrm{Re}(\widetilde {\gamma }^a) \lt 0\}\cup \{\widetilde {\gamma }^a : \textrm{Re}(\widetilde {\gamma }^a) = 0,\textrm{Im}(\widetilde {\gamma }^a) \gt 0\}, \end{equation}
a set that we use to construct the positive characteristic admittance. These correspond to purely forward-propagating waves, as well as forward-decaying waves that may propagate either forwards or backwards.
We must also now keep track of eigenvalue signs when defining eigenvectors. If we define the diagonal
$\widetilde {\boldsymbol{\varGamma }}^a$
matrix, and a general acoustic eigenvector
$(\widetilde {\boldsymbol{c}}_{\alpha ,1}^{a\pm },\widetilde {\boldsymbol{c}}_{\alpha ,2}^{a\pm })$
, there will exist matrices
$\widetilde { \unicode{x1D63E}}^{a\pm }$
, constructed from an ordering of either the ‘
$+$
’ lower eigenvectors or the ‘
$-$
’. The lower-block admittance equation is then
\begin{align} & \big(\widetilde { \unicode{x1D647}}_{\,3}^{\,a}\kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{\,a\pm} + \widetilde { \unicode{x1D647}}_1^a \big)\widetilde { \unicode{x1D63E}}^{a\pm } = \pm \widetilde { \unicode{x1D63E}}^{a\pm }\widetilde {\boldsymbol{\varGamma }}^{a},&&\textrm{so} && \kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{\,a\pm} = \left (\widetilde { \unicode{x1D647}}_{\,3}^{\,a}\right )^{-1}\left [-\widetilde { \unicode{x1D647}}_1^a \pm \widetilde { \unicode{x1D63E}}^{a\pm }\widetilde {\boldsymbol{\varGamma }}^{a}\big (\widetilde { \unicode{x1D63E}}^{a\pm }\big )^{-1}\right ], \end{align}
from which we deduce that the torsional entries
$\widetilde { \unicode{x1D647}}_1^a$
have caused a symmetry breaking between the positive and negative characteristic admittances. When it comes to computing the nonlinear admittances, we once again must take a transpose
\begin{equation} \kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{\,a\pm}\widetilde { \unicode{x1D647}}_{\,3}^{\,a} - \widetilde { \unicode{x1D647}}_1^a = \left (\widetilde { \unicode{x1D647}}_{\,3}^{\,a}\kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{\,a\pm} + \widetilde { \unicode{x1D647}}_1^a\right )^{\textrm{T}} = \pm \big (\widetilde { \unicode{x1D63E}}^{a\pm }\big )^{-{\textrm{T}}}\widetilde {\boldsymbol{\varGamma }}^{a}\big (\widetilde { \unicode{x1D63E}}^{a\pm }\big )^{\textrm{T}}, \end{equation}
meaning that the nonlinear characteristic admittances are given by
\begin{align} & \widetilde {\mathcal{Y}}_{\,\alpha \beta \gamma }^{\,ab\pm } = \pm \sum _{\zeta ,\delta ,\epsilon }\big (\widetilde { \unicode{x1D63E}}^{a\pm }\big )_{\zeta \alpha }^{-1}\big (\widetilde { \unicode{x1D63E}}^{(a - b)\pm }\big )_{\delta \beta }^{-1}\big (\widetilde { \unicode{x1D63E}}^{b\pm }\big )_{\epsilon \gamma }^{-1} \nonumber \\& \times \frac {\left [\big (\widetilde { \unicode{x1D63E}}^{a\pm } \big )^{{\textrm{T}}}\!\left [\widetilde {\mathcal{N}}_1^{\,\,ab} \big\langle \kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{(a - b)\pm },\kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{b\pm }\big\rangle + \widetilde {\mathcal{N}}_6^{\,\,ab} -\! \kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{\,a\pm}\widetilde {\mathcal{N}}_{\,\,7}^{\,\,ab}\big\langle \kern2.5pt \widetilde { \kern-2.5pt \unicode{x1D654}}^{(a - b)\pm }, \unicode{x1D644}\kern2pt\big\rangle \right ] \big \langle \widetilde { \unicode{x1D63E}}^{(a - b)\pm },\widetilde { \unicode{x1D63E}}^{b\pm }\big \rangle \right ]_{\zeta \delta \epsilon }}{\widetilde {\gamma }_\zeta ^a + \widetilde {\gamma }_\delta ^{a - b} + \widetilde {\gamma }_\epsilon ^b}. \end{align}
Appendix D. Viscosity in more detail
D.1. Numerical viscosity calculation
A more detailed derivation of the numerical viscosity term in (3.24) is presented here. Starting from (4.9), the summand may be split into partial fractions as
\begin{equation} \frac {a}{2(a - b)b} = \frac {1}{2(a - b)} + \frac {1}{2b} , \end{equation}
and when summed from
$-\infty$
to
$\infty$
this becomes
$-1/a$
. This can be folded into the other sum, so that the error becomes
\begin{equation} E^a = \frac {{\textrm{i}}\kern0.5pt\omega \beta_0 \sqrt {A_{{cs}}}M^2\textrm{e}^{{\textrm{i}}\kern0.5pt a\omega s}}{2(1 + \sigma )^2}\left (\sum _{\substack {b = a - {\textrm{sgn}}(a)a_{\textrm{max}},\\b \neq a}}^{{\textrm{sgn}}(a)a_{\textrm{max}}}\frac {1}{a - b} + \sum _{\substack {b = a - {\textrm{sgn}}(a)a_{\textrm{max}},\\b \neq 0}}^{{\textrm{sgn}}(a)a_{\textrm{max}}}\frac {1}{b}\right )\!. \end{equation}
The first sum is split as
\begin{equation} \sum _{\substack {b = a - {\textrm{sgn}}(a)a_{\textrm{max}},\\b \neq a}}^{{\textrm{sgn}}(a)a_{\textrm{max}}}\frac {1}{a - b} = \left (\sum _{b = a - {\textrm{sgn}}(a)a_{\textrm{max}}}^{2a - {\textrm{sgn}}(a)(a_{\textrm{max}} + 1)} + \sum _{\substack {2a - {\textrm{sgn}}(a)a_{\textrm{max}},\\ b \neq a}}^{{\textrm{sgn}}(a)a_{\textrm{max}}}\right )\frac {1}{a - b}, \end{equation}
with the second term vanishing due to symmetric limits and an odd summand about
$b = a$
. The second sum is treated similarly
\begin{equation} \sum _{\substack {b = a - {\textrm{sgn}}(a)a_{\textrm{max}},\\b \neq 0}}^{{\textrm{sgn}}(a)a_{\textrm{max}}}\frac {1}{b} = \left (\sum _{\substack {b = a - {\textrm{sgn}}(a)a_{\textrm{max}},\\b \neq 0}}^{{\textrm{sgn}}(a)a_{\textrm{max}} - a} + \sum _{{\textrm{sgn}}(a)(a_{\textrm{max}} + 1) - a}^{{\textrm{sgn}}(a)a_{\textrm{max}}}\right )\frac {1}{b}, \end{equation}
and this time the first term vanishes. Both of the remaining sums, with the correct substitution, become
\begin{equation} \sum _{b = a - {\textrm{sgn}}(a)a_{\textrm{max}}}^{2a - {\textrm{sgn}}(a)(a_{\textrm{max}} + 1)}\frac {1}{a - b} = \sum _{{\textrm{sgn}}(a)(a_{\textrm{max}} + 1) - a}^{{\textrm{sgn}}(a)a_{\textrm{max}}}\frac {1}{b} = \sum _{b = 0}^{|a| - 1}\frac {{\textrm{sgn}}(a)}{a_{\textrm{max}} - b}, \end{equation}
so then the error becomes
\begin{equation} E^a = \frac {{\textrm{i}}\kern0.5pt\omega \beta_0 \sqrt {A_{{cs}}}M^2\textrm{e}^{{\textrm{i}}\kern0.5pt a\omega s}}{(1 + \sigma )^2}\sum _{b = 0}^{|a| - 1}\frac {{\textrm{sgn}}(a)}{a_{\textrm{max}} - b} = \frac {|a|\omega \beta_0 MP_0^a}{1 + \sigma }\sum _{b = 0}^{|a| - 1}\frac {1}{a_{\textrm{max}} - b}. \end{equation}
The sum has upper and lower bounds in the form of integrals
\begin{equation} \int _0^{|a| - 1}\frac {{\rm d} x}{a_{\textrm{max}} - x} \lt \sum _{b = 0}^{|a| - 1}\frac {1}{a_{\textrm{max}} - b} \lt \int _0^{|a| - 1}\frac {{\rm d} x}{a_{\textrm{max}} - (x + 1)}. \end{equation}
We choose the lower bound as the approximation to the sum since the upper bound’s integrand is singular at its upper limit. Thus, evaluating the integral and including the viscous scale factor
$\nu _0$
, we arrive at (3.24).
D.2. Physical viscosity calculation
If we instead consider physical viscosity in one dimension, the (single) momentum equation picks up an extra term
\begin{equation} -{\textrm{i}}\kern0.5pt a\omega U^a + \frac {\partial P^a}{\partial s} - \left (\frac {4\mu }{3} + \zeta \right )\frac {\partial ^2U^a}{\partial s^2} = \frac {\partial Q^a}{\partial s}, \end{equation}
while the mass conservation equation is unchanged. We still wish to use the Riccati method, which requires our two equations to be first order in
$s$
; this can be achieved by substituting the
$O(M^2)$
expansion of
$\partial U^a/\partial s$
from the mass conservation equation here. That substitution results in
\begin{align} & -{\textrm{i}}\kern0.5pt a\omega U^a + \left (1 - {\textrm{i}}\kern0.5pt a\omega \left (\frac {4\mu }{3} + \zeta \right )\right )\frac {\partial P^a}{\partial s}\nonumber \\ & = \frac {\partial Q^a}{\partial s} - {\textrm{i}}\kern0.5pt a\omega \left (\frac {4\mu }{3} + \zeta \right )\frac {\partial }{\partial s}\left (\sum _b \beta_0 P^{a - b}P^b - Q^a\right )\!. \end{align}
The
$s$
-derivatives of
$O(M^2)$
terms may now have linear terms substituted into them. For notational ease, we define a modified frequency for each temporal index, and absorb both viscosities into a single constant
\begin{align} \varOmega ^a = \frac {{\textrm{i}}\kern0.5pt a\omega }{1 - {\textrm{i}}\kern0.5pt a\omega \tilde {\mu }}, && \tilde {\mu } = \frac {4\mu }{3} + \zeta , \end{align}
so then we have
\begin{equation} \frac {\partial P^a}{\partial s} - \varOmega ^aU^a = \frac {\partial Q^a}{\partial s} -\tilde {\mu } \beta_0\varOmega ^a\frac {\partial }{\partial s}\sum _bP^{a - b}P^b. \end{equation}
Defining the one-dimensional admittances in the usual way, we end up with the following two equations for them
\begin{align} \frac {{\rm d}\mathsf{Y}^a}{{\rm d}s} &= -\varOmega ^a(\mathsf{Y}^a)^2 + {\textrm{i}}\kern0.5pt a\omega , \\[-12pt] \nonumber \end{align}
\begin{align} \frac {{\rm d}\mathcal{Y}^{\,ab}}{{\rm d} s} &= -\mathcal{Y}^{\,ab}\big (\varOmega ^a \unicode{x1D654}^{\,a} + \varOmega ^{a - b} \unicode{x1D654}^{\,a - b} + \varOmega ^b \unicode{x1D654}^{\,b}\big )\notag \\&\quad -\frac {1}{\sqrt { \pi}R}\Bigg \{\big (\varOmega ^{a - b} \unicode{x1D654}^{\,a - b} + \varOmega ^b \unicode{x1D654}^{\,b}\big )\left (\frac { \unicode{x1D654}^{\,a}}{2} + \tilde {\mu } \beta_0\varOmega ^a\right )\notag \\&\quad + \frac {i\omega }{2}\big (a\gamma + a \unicode{x1D654}^{\,a - b} \unicode{x1D654}^{\,b} - b \unicode{x1D654}^{\,a} \unicode{x1D654}^{\,a - b} - (a - b) \unicode{x1D654}^{\,a} \unicode{x1D654}^{\,b}\big )\Bigg \}. \\[9pt] \nonumber \end{align}
Since we are working in one dimension, the admittances must be characteristic, and so we have the solutions
\begin{align} \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,a\pm } &= \pm \kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a} = \pm \exp \left (-\frac {{\textrm{i}}\kern0.5pt \pi}{4}\left [{\textrm{sgn}}(1 - {\textrm{i}}\kern0.5pt a\omega \tilde {\mu }) - 1\right ]\right )\sqrt {|1 - {\textrm{i}}\kern0.5pt a\omega \tilde {\mu }|}, \\[-12pt] \nonumber \end{align}
\begin{align} \overline {\mathcal{Y}}^{\,ab\pm } &= \mp \frac {1}{\sqrt { \pi}R}\Bigg \{\left (\varOmega ^{a - b}\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a - b} + \varOmega ^b\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,b}\right )\left (\frac {\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a}}{2} \pm \tilde {\mu } \beta_0\varOmega ^a\right )\nonumber\\ &\quad + \frac {{\textrm{i}}\kern0.5pt\omega }{2}\!\left (a\gamma + a\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a - b}\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,b}\! - b\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a}\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a - b}\!\! - (a - b)\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a}\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,b}\right )\!\Bigg \}\bigg / \nonumber \\ & \quad \left (\varOmega ^a\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a}\! + \varOmega ^{a - b}\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,a - b}\!\! + \varOmega ^b\kern1.5pt\overline { \kern-1.5pt \unicode{x1D654}}^{\,\,b}\right ). \\[9pt] \nonumber \end{align}
Assuming
$\tilde {\mu }$
is small, first-order expansions of these quantities may be found, although the nonlinear admittance only appears in
$O(M^2)$
terms to begin with, so its second term here will eventually be
$O(M^2\mu )$
, and thus discarded. Hence, we have
\begin{equation} \frac {{\rm d} P_0^a}{{\rm d} s} \approx {\textrm{i}}\kern0.5pt a\omega \left [P_0^a - \frac {{\textrm{i}}\kern0.5pt a\omega \beta_0}{2\sqrt { \pi} R}\sum _bP_0^{a - b}P_0^b\right ] -\frac {a^2\omega ^2\tilde {\mu }}{2}P_0^a, \end{equation}
suggesting that in order to stabilise a truncated sawtooth wave the physical viscosity would need to be around
$\tilde {\mu } \approx 2 \beta_0 M/a_{\textrm{max}}\omega$
.
Appendix E. Inverse exponential horn in more detail
This appendix contains a more detailed presentation of the analytical solution to the inverse exponential horn’s first antisymmetric mode equation, discussed earlier in § 4.5. By use of the substitution
$\sigma = \pi{\textrm{e}}^{2|m|s}/2|m|X_{\textrm{i}}$
, we reduced the problem to the modified Bessel equation
\begin{equation} \sigma ^2\frac {{\rm d}^2p}{{\rm d} \sigma ^2} + \sigma \frac {{\rm d} p}{{\rm d}\sigma } - (\nu ^2 + \sigma ^2)p = 0, \end{equation}
where
$\nu$
is defined by
$\nu ^2 = 1 - \omega ^2/4m^2$
, and takes purely real values or purely imaginary values, depending on whether
$\omega$
exceeds
$2|m|$
. If
$\nu$
is not a real integer, then the solution is given by
\begin{equation} p(\sigma ) = c_1{{I}}_\nu (\sigma ) + c_2{{K}}_{\nu }(\sigma ). \end{equation}
The velocity was eliminated in the derivation of (4.20), but can be retrieved from the fourth row of (4.19), since we have
\begin{align} \frac {{\rm d} p}{{\rm d} s} = {\textrm{i}}\kern0.5pt\omega \,u + \frac {X'}{X}p &&\implies &&u(\sigma ) = \frac {2|m|}{{\textrm{i}}\kern0.5pt\omega }\left [\sigma p'(\sigma ) + p(\sigma )\right ], \end{align}
and thus the velocity is given by
\begin{equation} u(\sigma ) = \frac {2|m|}{{\textrm{i}}\kern0.5pt\omega }\big[c_1\tilde {{{I}}}_\nu (\sigma ) + c_2\tilde {{{K}}}_{\nu }(\sigma )\big],\quad \text{where}\quad \tilde {{{I}}}_\nu (\sigma ) := \sigma {{I}}_\nu '(\sigma ) + {{I}}_\nu (\sigma ), \end{equation}
with
$\tilde {{{K}}}_\nu (\sigma )$
defined similarly. We then have an admittance given by
\begin{equation} Y(\sigma ) = \frac {u(\sigma )}{p(\sigma )} = \frac {2|m|}{{\textrm{i}}\kern0.5pt\omega }\frac {c_1\tilde {{{I}}}_\nu (\sigma ) + c_2\tilde {{{K}}}_{\nu }(\sigma )}{c_1{{I}}_\nu (\sigma ) + c_2{{K}}_{\nu }(\sigma )}. \end{equation}
Applying the boundary condition
$Y(\sigma _{{o}}) = \overline {Y}_{{o}}$
, we get
\begin{align} Y(\sigma ) = \dfrac {2|m|}{{\textrm{i}}\kern0.5pt\omega }\dfrac {\left [\tilde {{{K}}}_{\nu }(\sigma _{{o}}) - \dfrac {{\textrm{i}}\kern0.5pt\omega \overline {Y}_{{o}}}{2|m|}{{K}}_{\nu }(\sigma _{{o}})\right ]\tilde {{{I}}}_\nu (\sigma ) - \left [\tilde {{{I}}}_\nu (\sigma _{{o}}) - \dfrac {{\textrm{i}}\kern0.5pt\omega \overline {Y}_{{o}}}{2|m|}{{I}}_{\nu }(\sigma _{{o}})\right ]\tilde {{{K}}}_{\nu }(\sigma )}{\left [\tilde {{{K}}}_{\nu }(\sigma _{{o}}) - \dfrac {{\textrm{i}}\kern0.5pt\omega \overline {Y}_{{o}}}{2|m|}{{K}}_{\nu }(\sigma _{{o}})\right ]{{I}}_\nu (\sigma ) - \left [\tilde {{{I}}}_\nu (\sigma _{{o}}) - \dfrac {{\textrm{i}}\kern0.5pt\omega \overline {Y}_{{o}}}{2|m|}{{I}}_{\nu }(\sigma _{{o}})\right ]{{K}}_{\nu }(\sigma )}, \end{align}
where the characteristic admittance
$\overline {Y}_{{o}}$
is given by plugging the definition of
$X_{{o}}$
into (3.12)
\begin{equation} \overline {Y}_{{o}} = \overline {Y}(\sigma _{{o}}) = \exp \left \{-\frac {{\textrm{i}}\kern0.5pt \pi}{4}\left [{\textrm{sgn}}\left (1 - \frac {4m^2\sigma _{{o}}^2}{\omega ^2}\right ) -1 \right ]\right \}\sqrt {\left |1 - \frac {4m^2\sigma _{{o}}^2}{\omega ^2}\right |}. \end{equation}
Specifying the total pressure at the inlet, we match the overall remaining constant with (3.23) to get
\begin{align} p(\sigma ) = \dfrac {M\sqrt {X_{\textrm{i}}}}{2{{i}}}\dfrac {\left [\tilde {{{K}}}_{\nu }(\sigma _{{o}}) - \dfrac {{\textrm{i}}\kern0.5pt\omega \overline {Y}_{{o}}}{2|m|}{{K}}_{\nu }(\sigma _{{o}})\right ]{{I}}_\nu (\sigma ) - \left [\tilde {{{I}}}_\nu (\sigma _{{o}}) - \dfrac {{\textrm{i}}\kern0.5pt\omega \overline {Y}_{{o}}}{2|m|}{{I}}_{\nu }(\sigma _{{o}})\right ]{{K}}_{\nu }(\sigma )}{\left [\tilde {{{K}}}_{\nu }(\sigma _{{o}}) - \dfrac {{\textrm{i}}\kern0.5pt\omega \overline {Y}_{{o}}}{2|m|}{{K}}_{\nu }(\sigma _{{o}})\right ]{{I}}_\nu (\sigma _{\textrm{i}}) - \left [\tilde {{{I}}}_\nu (\sigma _{{o}}) - \dfrac {{\textrm{i}}\kern0.5pt\omega \overline {Y}_{{o}}}{2|m|}{{I}}_{\nu }(\sigma _{{o}})\right ]{{K}}_{\nu }(\sigma _{\textrm{i}})}. \end{align}
The boundary condition at the duct outlet
$\sigma _{{o}}$
can be rewritten as a vector orthogonality condition in two dimensions, i.e.
\begin{equation} \begin{pmatrix} c_1 \\ c_2 \end{pmatrix}\boldsymbol{\cdot }\left [\sigma _{{o}}\begin{pmatrix} {{I}}_\nu '(\sigma _{{o}}) \\ {{K}}_\nu '(\sigma _{{o}}) \end{pmatrix} + \left (1 - \frac {{\textrm{i}}\kern0.5pt\omega \overline {Y}_{{o}}}{2|m|}\right )\begin{pmatrix} {{I}}_\nu (\sigma _{{o}}) \\ {{K}}_\nu (\sigma _{{o}}) \end{pmatrix} \right ] = 0. \end{equation}
Note that
${{I}}_\nu (\sigma )$
and
${{I}}_\nu '(\sigma )$
are exponentially growing functions, while
${{K}}_\nu (\sigma )$
and
${{K}}_\nu '(\sigma )$
exponentially decay. Unscaled, these functions will therefore be of very different orders of magnitude at the outlet. This relation, however, effectively sets the vector
$(c_1,c_2)$
to be ‘asymptotically orthogonal’ to the Bessel function vector, meaning that the two terms in both the velocity and the pressure are instead of comparable order at the outlet. If we then move backwards from the outlet, the
${K}$
functions will grow rapidly while the
${I}$
functions decay, so for most of the duct the
${K}$
functions will dominate. For higher frequencies,
$\nu$
can become imaginary, causing
${K}$
to oscillate towards the inlet: these oscillations will pass through zero, resulting in pressure nodes, which cause singularities in the admittance. We note that for the inverse exponential horn, the lowest-mode admittance singularities correspond approximately to roots of
${{K}}_{{\textrm{i}}\sqrt {\omega ^2/4m^2 - 1}}( \pi{\textrm{e}}^{2|m|s}/2|m|X_{\textrm{i}})$
. Figure 17 shows these relationships in more detail, for a frequency low enough that oscillations do not yet occur.
Another feature of this geometry (or any geometry with width variation) is that turning points can occur, i.e. points where a mode can go from being cut-on to cutoff. These are significant because they correspond mathematically to the eigenvalue
$\overline {\gamma }_{\,\alpha}^{\,a}$
for that particular mode passing from the real axis to the imaginary axis, and therefore necessarily passing through zero. Since the splitting operators then involve the inverse of a matrix with a zero eigenvalue, this results in singularities in the forward and backward-going pressures at this point. By inverting (3.13) for
$\omega = \overline {\omega }_\alpha ^a, \lambda _\alpha = \alpha \pi$
, and
$X = X_{\textrm{i}}{\textrm{e}}^{2|m|s}$
, we note that the location of the turning point for mode
$(a,\alpha )$
is
\begin{equation} s_{\textit{tp}} = \frac {1}{2|m|}\log \left (\frac {a\omega X_{\textrm{i}}}{\alpha \pi}\right ). \end{equation}
Figure comparing the (real parts of the) logarithmically plotted modified Bessel functions (
${{I}}_\nu (\sigma )$
in blue,
${{K}}_\nu (\sigma )$
in red), along an inverse exponential horn of length
$4.5X_i$
and width decrease ratio of 4, for a frequency of 0.5/
$X_i$
. The r.m.s. pressure (according to this analytical solution) normalised by Mach number is logarithmically plotted in black.
Open access
